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A.  SCIENTIFIC  OBJECTIVES:  The  objective  of  this  research  program  is  to  apply 
model-based  signal  processing  methods  to  enhance  the  performance  of  computation 
electromagnetics  (CEM)  simulators  for  shipboard  antenna  design.  While  recent  advances 
in  CEM  algorithms  has  significantly  reduced  the  simulation  cost  of  modeling  complex 
radiation  and  scattering  phenomena,  real-world  engineering  design  and  optimization 
often  require  that  calculations  be  carried  out  repeatedly  over  large  parameter  spaces  such 
as  frequency  and  aspect  angle,  making  the  computation  cost  still  exceedingly  high.  The 
goal  of  this  research  is  to  apply  model-based  signal  processing  algorithms  to  CEM 
simulators  to  overcome  such  computational  bottleneck  and  achieve  higher  design 
throughput.  In  particular,  we  shall  address  the  issue  of  how  to  interpolate  and  extrapolate 
the  frequency  and  angular  behaviors  of  antenna  characteristics  based  on  a  sparse  set  of 
computed  data.  We  shall  develop  algorithms  to  extract  sparse  and  physical  models  of  the 
relevant  physics  embedded  in  CEM  data  to  facilitate  design  and  synthesis.  Finally,  we 
shall  explore  the  transportability  of  the  methodology  to  other  applications  such  as  radar 
signatures  and  wireless  channel  characteristics. 

B.  SUMMARY  OF  RESULTS  AND  SIGNIFICANT  ACCOMPLISHMENTS:  Our 

research  during  the  second  year  of  this  program  has  been  focused  on  three  topics.  First, 
we  continue  our  research  to  develop  an  effective  algorithm  to  extrapolate  the  frequency 
behavior  of  antenna  radiation  characteristics  on  a  complex  platform  from  a  limited  set  of 
computed  data.  A  frequency-dependent  time-of-arrival  model  is  utilized  to  more 


accurately  describe  the  scattering  physics  of  the  induced  current  on  the  platform.  We 
have  devised  and  implemented  a  robust  algorithm  to  estimate  the  model  parameters, 
including  the  additional  frequency-dependent  factors.  Second,  we  have  developed  an 
algorithm  based  on  the  adaptive  feature  extraction  approach  to  address  the  frequency 
interpolation  problem.  This  algorithm  has  also  been  extended  to  the  angular  domain  to 
achieve  two-dimensional  interpolation.  It  was  used  in  conjunction  with  the  fast  multipole 
code  FISC  to  predict  the  radar  image  of  the  benchmark  VFY-218  airplane  at  UHF  band 
with  excellent  results.  Third,  we  have  developed  a  methodology  to  parameterize  complex 
antenna  radiation  patterns  using  a  sparse  point  radiator  model.  Our  approach  is  based  on 
a  matching  pursuit  algorithm  and  the  concept  has  been  demonstrated  on  FISC-generated 
data  for  a  ship-like  platform.  The  detailed  descriptions  of  these  three  topics  are  described 
below.  They  are  followed  by  discussions  on  some  exploratory  efforts  related  to  this 
program. 

Model-based  frequency  extrapolation  of  antenna  radiation  patterns  on  complex 
platforms.  In  antenna  design  and  analysis,  the  mounting  platform  can  have  a  significant 
effect  on  the  antenna  radiation  characteristics.  However,  rigorous  solution  of  the  radiation 
problem  over  a  complex  platform  is  very  time  consuming,  and  the  computation 
complexity  increases  dramatically  as  the  frequency  increases.  During  the  first  year  of 
this  program,  we  developed  a  model-based  frequency  extrapolation  technique  with  which 
the  radiated  field  over  a  broad  band  of  frequencies  can  be  obtained  using  rigorously 
computed  results  at  low  frequencies.  Our  approach  entailed  fitting  the  currents  computed 
at  low  frequencies  to  a  time-of-arrival  model  and  determining  the  model  coefficients 
using  the  superresolution  algorithm  ESPRIT  [1].  The  currents  and  the  radiation 
characteristics  were  then  extrapolated  from  the  resulting  model.  Some  initial  results  were 
obtained  during  the  first  year  but  the  accuracy  of  the  extrapolated  results  was  not 
satisfactory. 

During  this  past  year,  we  have  significantly  improved  the  algorithm  by  adopting  an 
improved  frequency-dependent  model: 
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to  parameterize  the  induced  current  on  the  platform.  In  the  above  expression,  tn  is  the 
arrival  time  of  the  nth  incident  wave  and  an  and  a„  are  its  amplitude  and  frequency 

dependency,  respectively.  The  different  time-of-arrival  terms  correspond  to  the  different 
incident  wave  mechanisms  from  both  the  direct  antenna  radiation  and  higher-order 
scattering  from  other  parts  of  the  platform,  as  illustrated  in  Fig.  1.  These  different 
mechanisms  have  in  general  different  frequency  dependencies.  For  canonical  platform 
shapes,  their  exact  frequency  dependencies  are  known  through  the  geometrical  theory  of 
diffraction  (GTD)  [2].  However,  for  more  complex  structures,  they  must  be  determined 
numerically.  To  accurately  extrapolate  the  frequency  response  to  a  much  broader  range, 
the  accurate  estimation  of  the  frequency  factors  is  critical.  A  small  error  in  a  will  result  in 
dramatic  difference  in  amplitude  at  frequencies  in  the  extrapolated  region.  However,  the 
existing  superresolution  algorithms  based  on  eigenspace  decomposition  (e.g.  ESPRIT  and 
MUSIC)  cannot  be  directly  apply  to  this  more  realistic  model.  We  have  devised  a  pre¬ 
multiplication  scheme  in  conjunction  with  the  complex  time-of-arrival  estimation  from 
ESPRIT  to  determine  the  additional  frequency-dependent  factors  [22,24].  The 
performance  of  the  algorithm  in  the  presence  of  noise  has  been  evaluated  based  on 
simulated  data  and  errors  in  the  estimation  of  model  parameters  have  been  quantified. 
Our  results  show  that  the  method  is  quite  robust.  The  algorithm  has  been  applied  to 
extrapolate  the  induced  currents  and  radiation  patterns  in  both  2D  and  3D  antenna- 
platform  radiation  problems. 

As  an  example,  a  2D  structure  shown  in  Fig.  2(a)  is  considered.  The  platform  is  14m 
in  length  and  3m  in  height.  The  antenna  is  a  horizontal  line  source  placed  at  5m  above  the 
platform.  The  induced  current  on  the  platform  is  computed  from  0.1  to  0.5  GHz  at  21 
frequency  points  using  2D  method  of  moments  (MoM).  The  current  is  extrapolated  to  1 .3 
GHz  and  radiated  field  is  then  computed  based  on  the  extrapolated  current.  Both  the 
frequency-independent  and  frequency-dependent  time-of-arrival  models  are  used  to 
perform  the  extrapolation,  and  the  resulting  radiated  fields  at  the  observation  angle  of  0  = 
40°  are  plotted  in  Figs.  2(b)  and  2(c),  respectively.  Also  plotted  is  the  reference  MoM 
result  obtained  via  brute  force  computation.  The  primary  radiation  of  the  dipole  antenna 
is  not  included  in  the  plots  so  that  we  can  better  observe  the  secondary  radiation  from  the 
platform.  It  is  apparent  that  the  frequency  dependency  in  the  field  response  is  not 
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Fig.  1 .  Time-of-arrival  model  for  the  induced  current  at  point  P  accounts  for 

the  direct  incident  radiation  from  the  antenna  and  the  multiple  scattered 
waves  from  other  parts  of  the  platform. 
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Fig.  2(a).  2D  platform  geometry 


captured  by  the  frequency-independent  time-of-arrival  model,  while  the  field  predicted 
by  the  frequency  dependent  model  is  in  good  agreement  with  the  computed  result. 

Next,  we  look  at  a  3D  platform  shown  in  Fig.  3(a).  The  antenna  is  a  horizontal  dipole 
oriented  in  the  x  direction.  The  induced  current  is  computed  from  0.1  to  0.36  GHz  at  13 
frequencies  and  extrapolated  to  0.7GHz.  The  computation  is  carried  out  using  FISC,  a  3D 
MoM  code  based  on  the  fast  multipole  method  [3],  The  extrapolated  frequency  response 
at  the  observation  angle  </>e,  =  30°,  ^  =  -60°  is  plotted  in  Fig.  3(b).  Also  plotted  for 
comparison  is  the  reference  response  computed  by  FISC  via  brute  force.  The  major 
radiation  features  are  captured  by  the  extrapolation.  Figs.  3(c)  and  3(d)  show  the 
reference  and  extrapolated  radiation  patterns  as  functions  of  frequency  and  azimuth  angle 
when  the  elevation  angle  is  fixed  at  50°.  Good  qualitative  agreement  is  observed.  The 
correlation  index  between  the  two  figures  is  found  to  be  0.9980  in  the  sampled  region  and 
0.9742  in  the  extrapolated  region,  respectively.  The  computation  time  of  the  brute  force 
reference  results  is  about  50  hours  on  a  Pentiumll  400MHz  PC,  while  the  total 
computation  time  to  carry  out  the  electromagnetic  computation  in  the  sampled  region  and 
the  extrapolation  process  is  7  hours. 

Although  the  determination  of  model  parameters  is  more  complicated,  the  frequency 
dependent  model  show  significant  performance  improvement  over  the  frequency 
independent  model.  This  is  due  to  the  improved  modeling  of  the  scattering  physics.  Since 
the  current  required  for  the  extrapolation  is  only  computed  at  lower  frequencies,  large 
savings  in  computation  time  and  memory  can  be  achieved.  We  hope  to  further  extend 
this  methodology  to  model  the  near  fields  of  complex  platforms  for  problems  related  to 
antenna  coupling  and  radiation  hazard  evaluation  [17],  We  also  plan  to  investigate 
extrapolation  issues  related  to  non-perfect  conducting  structures. 

Model-based  frequency  interpolation  using  adaptive  feature  extraction.  We  have 
begun  to  investigate  the  interpolation  problem  for  CEM  data.  Although  the  potential 
payoff  of  an  interpolation  algorithm  is  not  as  great  as  extrapolation,  it  may  provide  a 
more  robust  way  to  achieve  computational  savings.  To  achieve  the  parameterization  in 
the  frequency  dimension,  we  again  utilize  the  multiple-arrival  model  for  the  induced 
current.  To  obtain  the  model  parameters  from  the  computed  data,  we  adopt  an  algorithm 
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Fig  3.  (b)  Frequency  response  predicted  by  the  frequency  dependent 
model  at  (bei  =  30°,  <baz-  -60°. 


termed  adaptive  feature  extraction  (AFE).  AFE  can  be  considered  as  a  generalization  of 
the  CLEAN  algorithm  [4],  and  is  similar  to  the  adaptive  joint  time-frequency  technique 
[5-7]  and  the  matching  pursuit  algorithm  [8]  in  the  iterative  manner  it  performs  the 
parameterization.  We  have  previously  applied  it  to  construct  the  inverse  synthetic 
aperture  radar  (ISAR)  image  from  radar  measurement  data  that  was  undersampled  in  the 
aspect  dimension  [18].  The  essential  idea  of  the  AFE  algorithm  is  to  search  and  extract 
out  individual  scattering  features  from  the  data  set  one  at  a  time.  During  each  iteration 
the  strongest  feature  is  identified  and  removed  from  the  original  data.  The  procedure  is 
then  iterated  until  the  data  is  well  parameterized  by  the  feature  set.  In  this  manner,  the 
interference  between  different  scattering  features,  which  is  significant  for  undersampled 
data,  can  be  largely  avoided.  Since  the  features  in  our  model  are  exponential  functions  of 
frequency  and  aspect,  we  use  random  sampling  during  the  collection  of  the  original 
computation  data  in  order  to  avoid  the  ambiguity  in  selecting  the  strongest  feature.  This 
approach  is  similar  to  the  random  array  concept  that  uses  highly  thinned  but  randomly 
spaced  elements  to  avoid  grating  lobes  [9]. 

We  have  also  extended  this  algorithm  to  interpolate  2-dimensional  frequency  and 
angle  data  [20,  26].  In  this  case,  a  more  generalized  time-of-arrival  and  angle-of-arrival 
model  is  used  to  parameterize  the  current.  We  have  applied  the  2D  algorithm  to  predict 
the  ISAR  image  of  the  benchmark  VFY218  airplane  [10].  The  fuselage  length  of  the 
airplane  is  15.33  meters  and  the  maximum  width  measured  from  the  wing  tips  is  8.90 
meters.  To  generate  its  ISAR  image  at  a  center  frequency  of  400  MHz  with  bandwidth 
267  MHz,  we  use  the  fast  solver  FISC  on  a  Pentium  II  450MHz  computer.  The  total 
number  of  sampling  points  must  be  at  least  40  frequency  points  by  40  aspect  points 
within  the  40-degree  aperture  to  satisfy  the  Nyquist  sampling  criterion.  The  resulting 
range  and  cross  range  resolution  is  about  half  a  meter.  Since  the  calculation  for  one 
single  frequency-aspect  point  takes  about  3  hours  (with  about  80,000  unknowns  in  the 
moment  equation),  the  total  computation  time  would  be  200  days  if  we  use  brute-force 
calculation  to  generate  the  data.  Based  on  a  5:1  undersampling  rate,  we  only  compute  the 
scattering  problem  at  62  randomly  sampled  points  and  use  the  AFE  interpolation  scheme 
to  interpolate  the  data  to  all  40x40=1600  points.  The  computation  time  is  about  8  days. 
The  resulting  ISAR  image  at  the  130-degree  (from  nose-on)  look  angle  is  plotted  in  Fig. 


4(a).  For  comparison,  Fig.  4(b)  shows  the  ISAR  image  constructed  from  the  chamber 
measurement  data  for  the  same  look  angle.  The  target  outline  is  overlaid  on  the 
measurement  image.  The  measurement  was  carried  out  at  the  US  Navy  China  Lake 
facility  on  a  1:30  scaled  model  at  the  frequency  band  of  8  to  16  GHz.  Comparing  Fig. 
4(a)  to  Fig.  4(b),  we  find  that  all  the  features  in  the  measurement  image  are  well  predicted 
in  the  simulated  image  from  using  FISC  and  interpolation. 

The  AFE  approach  has  been  found  to  be  stable  and  robust.  Sufficient  accuracy  in  the 
predicted  image  features  can  be  achieved  even  when  the  original  computed  data  is 
sampled  at  5:1  below  the  Nyquist  criterion  in  either  frequency  or  aspect.  Therefore,  the 
expected  time  savings  in  using  this  approach  is  about  25:1  in  comparison  to  the  brute- 
force  computation.  The  AFE  algorithm  does  involve  exhaustive  search  and  must  be 
carried  out  for  every  current  element  on  the  target.  However,  the  time  consumed  in  the 
interpolation  is  still  relatively  insignificant  when  compared  to  the  electromagnetic 
computation  time.  This  error  sources  that  limit  the  dynamic  range  of  the  final  interpolated 
data  are  currently  being  investigated.  The  approach  will  also  be  further  explored  for 
application  in  the  antenna  placement  problem. 

Extraction  of  sparse  representation  of  antenna  radiation  data.  The  previous 
topics  are  focused  on  the  utilization  of  model-based  processing  to  extrapolate  or 
interpolate  CEM  data  for  the  “forward”  solution  of  electromagnetic  radiation  and 
scattering  problems  involving  complex  platforms.  An  equally  important  problem  from 
the  design  perspective  is  the  “inverse”  algorithm  of  spatially  pinpointing  the  locations  of 
platform  scattering  based  on  CEM  data.  Toward  this  end,  we  had  previously  developed  a 
method  to  extract  sparse  model  of  the  antenna  radiation  pattern  on  a  complex  platform 
[11], [12].  This  representation  is  based  on  a  point  radiator  model  that  describes  the 
radiation  pattern  by  a  collection  of  “radiation  centers”  on  the  platform.  The  methodology 
for  obtaining  the  radiation  center  model  entails  first  generating  the  3-D  antenna  synthetic 
aperture  radar  (ASAR)  imagery  of  the  platform,  and  then  parameterizing  the  resulting 
image  by  a  collection  of  point  radiators  via  the  CLEAN  algorithm  [27].  Once  such  a 
representation  is  obtained,  we  can  rapidly  reconstruct  antenna  radiation  patterns  over 
frequencies  and  aspects  with  good  fidelity.  Thus  such  a  model  can  be  used  for  real-time 


0 


-5  0  5  10  15 


Range 


Fig.4(a)  ISAR  image  of  VFY-218  at  130  degrees  from  nose-on 
generated  from  interpolated  result  using  AFE  with  62 
FISC-computed  points. 
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Fig.4(b)  ISAR  image  of  VFY-218  at  130  degrees  from  nose-on 
generated  from  chamber  measurement  data. 


reconstruction  of  complex  antenna  patterns  in  high-level  system  simulations. 
Furthermore,  the  resulting  radiation  center  information  can  be  used  to  pinpoint  cause- 
and-effect  in  platform  scattering  and  provide  important  design  guidelines  for  antenna 
placement  and  optimization.  However,  the  method  was  based  on  a  Fourier-based 
algorithm  that  relied  on  a  number  of  small-angle,  small-bandwidth  approximations. 
Furthermore,  the  concept  was  only  demonstrated  using  high-frequency  ray-tracing 
simulations  and  not  more  rigorous  CEM  data. 

During  the  past  year,  we  have  overcome  the  above  deficiencies  by  developing  a  more 
generalized  matching  pursuit  algorithm  in  the  frequency-aspect  domain  based  on  a  point 
radiator  basis  [25],  We  have  also  demonstrated  the  feasibility  of  extracting  a  sparse 
model  of  the  antenna-platform  interaction  using  CEM  data  from  FISC.  The  matching 
pursuit  algorithm  is  implemented  based  on  the  following  radiation  center  basis: 


Es(k  $  (j))  =  Ae-jkroejk(xosin6cos<t>+yosinGsintt>+zoCos0)  _  Ae~jkr0  ej(kxx0+kyy0+kzzQ) 


(2) 

where  k  is  the  wave  number,  (xo,  yo,  zo)  is  the  location  of  the  radiation  center.  The  origin 
of  the  above  basis  is  illustrated  in  Fig.  5.  Note  that  the  phase  factor  in  (2)  accounts  for 
the  path  delay  from  the  antenna  to  the  point  scatterer  on  the  platform,  and  then  to  the 
observation  direction  in  the  far  field.  To  speed  up  the  parameterization  time  of  the 
matching  pursuit  algorithm,  we  estimate  the  location  of  the  radiation  centers  by  utilizing 
the  Fourier-based  ASAR  algorithm.  The  point  with  the  highest  intensity  is  first  located  in 
the  ASAR  image  and  its  amplitude  and  coordinates  are  determined  to  serve  as  an  estimate 
of  the  strongest  radiation  center.  We  next  zoom  in  on  the  precise  location  of  the  radiation 
center  via  a  fine  search.  We  then  subtract  the  contribution  of  this  radiation  center  from 
the  total  radiated  field  and  iterate  the  process  until  the  energy  in  the  residual  signal  has 
reached  a  sufficiently  small  level. 

As  an  example,  the  radiation  center  model  is  extracted  from  the  computed  radiation 
data  from  FISC  for  the  ship-like  structure  shown  earlier  in  Fig.  3(a).  The  radiated  field  is 
computed  at  a  center  frequency  of  1 .0  GHz  and  a  bandwidth  of  500  MHz.  The  center 
observation  angle  is  0e/  =  40°,  0^  =  50°  and  the  angular  range  is  about  23°  in  both  azimuth 
and  elevation  angles.  The  first  20  extracted  radiation  centers  are  plotted  in  Fig.  6(a).  The 
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Fig.  6(a).  First  20  radiation  centers  extracted  using  the  matching 
pursuit  algorithm 
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Fig.  6(b)  Original  radiated  field  computed  by  FISC  over  a 
bandwidth  of  500  MHz  and  23°  azimuth  angle. 
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Fig.  6(c)  The  radiated  field  reconstructed  from  the  first  20 
extracted  radiation  centers 


strength  of  the  radiation  centers  is  represented  by  different  colors.  We  observe  that  the 
dominant  platform  scattering  comes  from  the  edge  diffraction  from  the  bow  of  the  ship 
platform  and  the  comer  structure  formed  by  the  cylinder  and  plate.  Note  that  a  few 
radiation  centers  due  to  the  edge  diffraction  are  slightly  off  the  platform.  This  is  due  to 
the  limited  resolution  of  the  matching  pursuit  algorithm.  Once  the  sparse  representation 
is  generated,  the  radiated  field  can  be  easily  reconstructed  using  the  radiation  centers.  The 
original  radiated  field  at  kz  -  0  is  plotted  in  Fig.  6(b)  as  a  function  of  the  kx  and  ky,  where 
the  wave  numbers  kx,  ky  and  kz  are  defined  in  (2).  The  field  reconstructed  from  the  first  20 
radiation  centers  is  shown  in  Fig.  6(c).  The  two  patterns  match  well  over  both  frequency 
and  angle.  The  correlation  index  between  the  two  figures  is  found  to  be  0.958. 

We  have  demonstrated  that  radiation  patterns  from  complex  platforms  can  be  well 
represented  by  a  very  sparse  set  of  radiation  centers.  The  resulting  model  can  be  used  to 
spatially  pinpoint  cause-and-effect  on  the  platform  and  could  play  an  important  role  in 
antenna  placement  and  optimization.  We  shall  continue  to  refine  this  algorithm  and 
examine  more  complex  topside  structures.  We  also  plan  to  utilize  this  model  for 
addressing  the  antenna  coupling  issue  once  the  CEM  capability  for  simulating  antenna 
coupling  characteristics  is  in  place. 

Other  exploratory  topics.  Several  additional  studies  have  been  carried  out  during 
the  past  year  and  are  described  below.  They  serve  as  possible  launching  points  into  more 
relevant  and  focused  efforts  for  the  present  program  in  the  next  year.  First,  we  have 
begun  to  investigate  the  resonant  behaviors  of  antenna-platform  configurations.  This 
topic  is  motivated  by  the  interest  of  the  Navy  SPAWAR  Center  in  HF  antennas  where 
ship  body  resonances  can  dominate  the  antenna  frequency  characteristics.  In  the  lower 
frequency  regime,  the  scattering  phenomenology  differs  from  ray-optical  characteristics. 
Consequently,  the  time-of-arrival  model  we  have  utilized  successfully  may  no  longer  be 
efficient  in  this  regime.  We  have  carried  out  a  preliminary  study  to  understand  the 
scattering  phenomenology  in  these  configurations.  Furthermore,  we  are  exploring  a 
hybrid  scheme  to  achieve  a  sparse  parameterization  of  the  data  using  a  combination  of  (i) 
a  rational  function  model  for  those  regions  on  the  structure  dominated  by  resonant 
effects,  and  (ii)  a  exponential  model  for  those  regions  on  the  structure  dominated  by  ray 


optical  phenomena.  We  hope  to  devise  a  robust  parameterization  procedure  so  that  both 
types  of  physical  mechanisms  can  be  accurately  described. 

An  important  issue  related  to  the  study  of  resonant  region  using  an  iterative  CEM 
solver  like  FISC  is  that  when  the  platform  exhibits  strong  resonance  effect,  the  iteration 
number  required  for  an  accurate  solution  can  become  large.  A  good  preconditioner  for 
the  moment  matrix  is  needed  to  alleviate  this  problem.  We  have  been  investigating  the 
use  of  wavelet  packet  basis  for  the  sparsification  of  moment  matrix  [13,14].  We  have 
recently  extended  this  work  by  devising  an  approximate-inverse  preconditioner  to  the 
moment  equation  using  the  wavelet  representation  [15,23].  We  plan  to  further  investigate 
this  topic,  which  may  become  important  for  simulating  near-resonant  structures. 

C.  FOLLOW-UP  STATEMENT: 

In  the  coming  year,  we  will  continue  our  research  along  the  three  research  topics 
outlined  above  by:  (i)  extending  our  work  on  frequency  extrapolation  to  near-field 
radiation  characteristics,  (ii)  extending  the  frequency-aspect  interpolation  algorithm  for 
multiple  frequency-antenna  position  study,  and  (iii)  refining  the  radiation  center 
extraction  algorithm.  In  addition,  we  plan  to  initiate  several  new  research  topics  in:  (i) 
applying  the  new  FISC  code  for  mixed  wire-facet  bases  for  antenna  self-impedance 
study,  (ii)  investigating  antenna  mutual  coupling  issues  and  exploring  the  usage  of  CEM 
data  in  array  calibration,  and  (iii)  devising  an  improved  model  for  describing  both  ray- 
optical  and  resonance  phenomena  in  antenna-platform  interactions. 
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ABSTRACT:  An  angular  extrapolation  technique  is  proposed  to  acceler¬ 
ate  the  multiangle  RCS  calculation  of  iterative  MoM  solvers.  Based  on 
high-frequency  electromagnetic  phenomenology ,  an  exponential  model  is 
proposed  to  parameterize  the  induced  currents  on  the  target.  The  model 
parameters  are  extracted  from  MoM  results  at  a  limited  number  of 
angular  samples  using  the  superresolution  algorithm  ESPRIT.  The  RCSs 
at  a  wide  extent  of  angles  are  then  extrapolated.  This  approach  is  tested 
on  two  numerical  examples.  Comparison  with  the  bistatic  approximation 
shows  the  improved  accuracy  of  the  algorithm.  ©  1999  John  Wiley  & 
Sons,  Inc.  Microwave  Opt  Technol  Lett  20:  229-233,  1999. 

Key  words:  angular  extrapolation ;  method  of  moments;  superresolution 

I.  INTRODUCTION 

It  is  well  known  that  when  the  method  of  moments  (MoM) 
equation  is  solved  using  direct  matrix  inversion  or  LU  decom¬ 
position,  multiple  right-hand-side  (RHS)  solutions  can  be 
generated  efficiently.  Unfortunately,  the  computation  com¬ 
plexity  of  OiN3)  for  such  direct  solvers  is  usually  too  high 
when  the  target  size  is  bigger  than  a  few  wavelengths.  To 
overcome  this  difficulty,  many  iterative  MoM  solvers  have 
been  developed  by  using  iterative  matrix-vector  multiplica¬ 
tion.  Among  them,  the  fast  multipole  method  (FMM)  [1]  and 
the  multilevel  fast  multipole  method  (MLFMM)  [2]  can  re¬ 
duce  the  computational  complexity  of  the  matrix-vector  mul¬ 
tiply  operation  from  0(N2)  to  0(NL5)  and  OiN  log  A/), 
respectively.  However,  for  applications  such  as  generating 
multiangle  radar  cross  section  (RCS)  data,  these  iterative 
procedures  have  to  be  repeated  for  each  RHS  or  look  angle. 
There  are  some  existing  techniques  to  speed  up  iterative 
matrix  solvers  for  multiple  RHS  problems  [3,  4],  Nevertheless, 
the  computational  gain  usually  comes  at  the  expense  of  more 
memory  requirement.  In  this  paper,  we  propose  an  angular 
extrapolation  scheme  to  extrapolate  angular  RCS  data  from 
the  solutions  at  a  few  densely  sampled  angles.  These  solu¬ 
tions  can  be  obtained  by  iterative  MoM  solvers,  either  with  or 
without  the  accelerated  multiple  RHS  algorithms. 

An  often-used  angular  extrapolation  (or  interpolation) 
approach  for  RCS  is  the  so-called  bistatic  approximation 
[5-7],  which  is  based  on  the  physical  optics  assumption. 
Under  this  assumption,  the  monostatic  scattered  data  are 
approximated  by  the  bistatic  scattered  data  at  nearby  angles. 


MICROWAVE  AND  OPTICAL  TECHNOLOGY  LETTERS  /  Vol.  20,  No.  4,  February  20  1999 


Because  bistatic  data  are  relatively  inexpensive  to  compute, 
additional  monostatic  data  can  be  readily  generated  using 
this  approach.  However,  since  physical  optics  does  not  take 
into  account  multiple  interactions  among  different  parts  of 
the  target,  the  bistatic  approximation  breaks  down  when 
applied  to  scattering  data  of  targets  that  are  dominated  by 
multiple  scattering  mechanisms. 

The  angular  extrapolation  algorithm  we  will  present  here 
is  based  on  a  multiple  scattering  model.  We  parameterize  the 
current  on  the  target,  which  is  available  from  the  MoM 
solution,  using  an  exponential,  multiple-excitation  model.  The 
parameters  of  this  model  are  extracted  from  a  few  angular 
samples  via  the  ESPRIT  superresolution  algorithm  [8],  Once 
these  parameters  are  obtained,  the  angular-dependent  model 
of  the  induced  current  at  each  point  on  the  target  is  con¬ 
structed.  Thus,  the  induced  current  over  a  wider  angular 
extent  can  be  extrapolated  and  the  multiangle  RCS  can  be 
calculated.  This  approach  can  be  considered  as  a  dual  of  the 
model-based  frequency  extrapolation  that  we  have  devised 
earlier  [9].  In  the  following  section,  the  multiple-excitation 
model  will  be  described  in  detail.  In  Section  III,  numerical 
examples  are  presented  to  illustrate  the  performance  of  our 
extrapolation  algorithm,  and  comparisons  are  made  against 
the  bistatic  approximation. 

II.  MODEL-BASED  EXTRAPOLATION 

The  key  to  the  success  of  any  model-based  extrapolation  is  a 
good  model  that  coincides  with  the  underlying  physical  mech¬ 
anisms  [9,  10].  We  postulate  that  in  high-frequency  scattering, 
the  induced  current  at  each  point  S  on  the  target  is  excited 
by  multiple-arriving  scattered  waves,  as  shown  in  Figure  1.  If 
we  denote  the  down-range  direction  with  respect  to  the 
incident  wave  as  x  and  the  cross-range  direction  as  y,  the 
current  at  S  can  be  written  as 

N 

7(5)=  d,-x,  +  l,  (1) 

/=  I 

where  N  is  the  number  of  scattered  waves  arriving  at  S  and 
k  is  the  free-space  wavenumber.  In  the  above  definition  of 
the  path  length  dh  we  let  y,)  be  the  first  hit  point  on  the 
target  due  to  the  incident  wave,  and  let  /,  be  the  total 
intermediate  path  length  of  the  multiple-scattering  mecha¬ 
nism.  Note  that  the  two  mechanisms  illustrated  in  Figure  1 
are  both  two-bounce  mechanisms.  For  a  single-bounce  physi¬ 
cal-optics  mechanism,  /,  =  0  and  (jc,,  yy)  simply  corresponds 
to  the  point  S . 


Figure  1  Multiple-arrival  model  of  the  scattering  mechanisms 


We  now  make  a  key  assumption  about  the  current  model 
as  the  incident  angle  is  varied.  We  assume  that  all  intermedi¬ 
ate  scattering  points  for  each  mechanism  remain  fixed  as  the 
incident  angle  is  varied,  as  illustrated  in  Figure  1.  This 
assumption  has  been  found  to  be  fairly  accurate  for  ray-opti¬ 
cal  fields  under  small  angular  variation  [11],  and  leads  us  to 
the  following  current  model  as  a  function  of  incident  angle: 

N 

7(0, S)  =  ^Ai(9)e-2jk^cose+y‘sin&+l'\  (2) 

i=  I 

If  we  further  use  the  small-angle  approximation  cos  0=1 

and  sin  0=0,  and  assume  that  is  independent  of  angle, 

then 

N 

/(0,S)=  (3a) 

/=  i 
N 

=  (3b) 

/=1 

We  observe  that  the  above  equation  is  in  the  form  of  a 
sum-of-exponential  model,  with  linear  phase  dependence  with 
respect  to  the  incident  angle. 

With  the  above  multiple-excitation  model  in  hand,  we  set 
out  to  extract  the  parameters  (^',y;)  at  each  point  on  the 
target  from  a  limited  number  of  angular  samples  calculated 
using  the  MoM.  We  apply  the  superresolution  algorithm 
ESPRIT  [9]  for  this  purpose  due  to  its  robustness  in  the 
presence  of  noise.  ESPRIT  is  based  on  the  data  model 

N 

F{8J=  £  +  n{ 0J,  9m  =  6„  d2, . . . ,  0M 

i=  1 

(4) 

where  n(-)  denotes  additive  white  Gaussian  noise  and  the 
samples  are  uniformly  spaced.  If  the  data  sequence  obeys  this 
ideal  model  exactly  and  the  number  of  sampling  points  M  is 
infinite,  ESPRIT  can  estimate  N  and  resolve  each  A)  and  y, 
without  any  error.  For  finite-length  data,  the  minimum  num¬ 
ber  of  samples  to  perform  the  estimation  is  M  >  2N  +  1, 
and  the  accuracy  of  the  estimated  parameters  will  depend  on 
the  length  of  the  available  data.  At  frequencies  above  target 
resonance,  the  first  few  interactions  dominate  the  scattering 
phenomenon,  and  only  a  small  N  is  needed  to  adequately 
model  the  observed  data.  Therefore,  by  running  ESPRIT 
using  the  MoM-generated  current  at  M  angles,  the  angular- 
dependent  current  model  in  (3b)  can  be  obtained.  Note  that 
the  ESPRIT  processing  is  carried  out  at  each  point  on  the 
target.  Once  such  a  model  is  found,  we  can  extrapolate  the 
data  to  determine  the  current  at  other  angles.  The  total 
scattered  field  versus  angle  can  thus  be  calculated  by  integrat¬ 
ing  the  extrapolated  current.  For  large  angle  extrapolation, 
we  find  that  it  is  more  accurate  to  use  the  original  sin  0 
factor  in  place  of  the  0  term  in  the  phase  of  the  current 
model  in  (3b).  To  both  satisfy  the  Nyquist  sampling  criterion 
for  ESPRIT  and  achieve  the  maximum  calculation  efficiency, 
the  angular  sampling  interval  is  normally  chosen  as  7r/2Aymax, 
where  ymax  is  the  maximum  target  dimension  in  cross  range. 

III.  NUMERICAL  RESULTS 

To  demonstrate  the  performance  of  the  angular  extrapolation 
algorithm,  two  numerical  examples  are  presented  based  on 
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Angle  (degrees) 

Figure  2  Comparison  between  the  single-point  bistatic  approximation  and  the  ESPRIT  extrapolation  results  for  the  dihedral 


Figure  3  Comparison  between  the  eight-point  bistatic  interpolation  results  and  the  ESPRIT  extrapolation  results  for  the  dihedral 
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2-D  MoM  results.  The  first  example  is  a  2-D  dihedral  corner 
reflector  (shown  in  the  inset  of  Fig.  2)  that  is  40  cm  along  the 
horizontal  face  and  60  cm  along  the  vertical  face.  The  fre¬ 
quency  of  the  incident  wave  is  3  GHz.  The  RCS  data  versus 
angle  are  calculated  by  the  MoM  code  between  25  and  65° 
(measured  from  the  horizontal  face).  The  data  are  plotted  as 
the  solid  curve  in  Figure  2.  A  relatively  isotropic  scattering 
pattern  is  observed  within  this  angular  range.  Using  eight 
points  evenly  sampled  between  43  and  50°.  we  carry  out  our 
extrapolation  algorithm.  The  model  order  N  used  is  3.  The 
extrapolated  RCS  data  are  plotted  as  the  dashed  curve  in 
Figure  2.  The  ESPRIT  result  coincides  well  with  the  exact 
result,  unless  the  extrapolation  angles  are  too  far  removed.  If 
we  use  the  bistatic  approximation  with  the  incident  angle  at 
45°  to  approximate  the  monostatic  RCS  in  this  range,  the 
RCS  data  fall  apart  rapidly  away  from  the  center  angle,  as 
shown  by  the  dash-dotted  line  in  Figure  2.  This  is  because 
the  dominant  double-bounce  mechanism  is  not  included  in 
the  bistatic  model.  Even  if  we  calculate  the  bistatic  field  at 
eight  points  evenly  distributed  in  the  angular  aperture  and 
linearly  interpolate  the  monostatic  RCS  in  between,  the 
result  still  deviates  from  the  exact  result  once  the  interpola¬ 
tion  angle  is  away  from  those  calculated  angles.  This  is 
illustrated  in  Figure  3. 

A  second  example  is  a  (circular  cylinder)-(plate)  scatterer 
shown  in  Figure  4.  The  diameter  of  the  cylinder  is  42  cm, 
and  the  length  of  the  plate  is  2  m.  The  distance  between 
the  origin  of  the  cylinder  and  the  plate  is  62  cm.  The  calcula¬ 
tion  frequency  is  3  GHz.  In  this  structure,  strong  multiple¬ 
scattering  mechanisms  between  the  cylinder  and  the  plate 
dominate  the  backscattering.  The  exact  RCS  data  calculated 
by  the  MoM  are  plotted  as  the  solid  curve  in  both  Figures  5 
and  6.  They  consist  of  a  total  of  81  points  from  25  to  65°. 
Only  eight  points  between  43.5  and  47°  are  used  for  the 
ESPRIT  extrapolation.  The  model  order  N  used  is  3.  The 
result  is  plotted  as  the  dashed  curve  in  Figure  5,  which  agrees 


42  cm 


Figure  4  Geometry  of  the  (circular  cylinder)-(plate)  target 

well  with  the  exact  result.  It  is  observed  that  some  deviations 
occur  when  the  extrapolation  angle  is  beyond  the  range 
35-55°.  Using  eight  points  evenly  distributed  in  the  angular 
aperture  [25°,  65°],  the  bistatically  interpolated  result  is  plot¬ 
ted  as  the  dashed  curve  in  Figure  6.  It  can  be  seen  that  the 
result  deviates  significantly  from  the  exact  result  between  any 
two  calculated  angles.  The  beatings  in  the  RCS  curve  are  not 
well  predicted. 

IV.  CONCLUSION 

A  model-based  extrapolation  algorithm  has  been  proposed  to 
generate  multiangle  RCS  for  iterative  MoM  solvers.  An  expo¬ 
nential  model  is  used  to  incorporate  the  multiple-scattering 
mechanisms.  The  parameters  of  the  exponential  model  are 
extracted  using  the  ESPRIT  superresolution  algorithm.  Once 
the  model  is  found,  it  can  be  extended  to  predict  the  RCS  at 
a  wider  range  of  angles.  Numerical  examples  show  that  the 
extrapolation  technique  has  robust  performance  for  targets 
dominated  by  multiple-scattering  mechanisms,  while  the  stan¬ 
dard  bistatic  approximation  usually  gives  poor  extrapolation/ 
interpolation  results  for  this  type  of  target.  The  accuracy  of 


Angle  (degrees) 

Figure  5  Comparison  between  the  ESPRIT  extrapolation  result  and  the  exact  result  for  the  cylinder  plate 
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Figure  6  Comparison  between  the  bistatic  approximation  result  and  the  exact  result  for  the  cylinder  plate 


the  present  extrapolation  algorithm  is  limited  by  the  validity 
of  the  model,  in  which  the  small-angle  approximation  is 
implicitly  assumed.  It  is  expected  that  multiplicative  speedup 
in  time  can  be  achieved  when  this  algorithm  is  combined  with 
the  existing  acceleration  schemes  for  multiple  RHSs  [3,  4]. 
Finally,  the  angular  extrapolation  algorithm  can  be  combined 
with  previously  proposed  frequency  extrapolation  to  generate 
2-D  frequency-aspect  data  for  radar  image  formation.  These 
topics  are  currently  under  study. 
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ACSAR- ANTENNA  COUPLING  SYNTHETIC  APERTURE 
RADAR  IMAGING  ALGORITHM 

C.  Ozdemir  and  H.  Ling 
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The  University  of  Texas  at  Austin 
Austin,  TX  78712-1084,  USA 

Abstract-A  synthetic  aperture  radar  imaging  technique  called  AC¬ 
SAR  for  antenna  coupling  scenarios  is  introduced.  It  is  shown  that  an 
ACSAR  image  of  a  platform  can  be  formed  by  inverse-Fourier  trans¬ 
forming  the  multi-frequency,  multi-spatial  coupling  data  between  two 
antennas.  Furthermore,  we  present  a  fast  ACSAR  imaging  algorithm 
that  is  specifically  tailored  to  the  shooting  and  bouncing  ray  (SBR) 
technique.  The  fast  algorithm  is  shown  to  reduce  the  total  simulation 
time  by  several  orders  of  magnitude  without  significant  loss  of  fidelity. 
Finally,  a  sparse  representation  of  ACSAR  imagery  is  introduced  by 
extracting  the  point  radiators  in  the  image.  By  parameterizing  the 
ACSAR  image,  it  is  possible  to  reconstruct  the  3-D  ACSAR  image 
and  the  3-D  frequency-spatial  data  with  a  very  sparse  set  of  radiation 
centers. 


1.  INTRODUCTION 

The  electromagnetic  coupling  between  antennas  on  a  complex  platform 
is  an  important  issue  in  antenna  design.  For  an  antenna  designer,  it 
is  useful  to  identify  the  areas  on  the  platform  that  cause  strong  cou¬ 
pling  between  the  antennas.  In  this  work,  we  set  out  to  develop  an 
imaging  algorithm  to  pinpoint  the  dominant  scattering  locations  on 
the  platform  that  give  rise  to  antenna  interactions  from  the  coupling 
data  between  antennas.  Our  approach  to  this  problem  is  based  on  the 
Inverse  Synthetic  Aperture  Radar  (ISAR)  concept.  ISAR  imaging  is  a 
standard  technique  to  map  the  locations  of  dominant  scattering  off  a 
target  based  on  the  multi-frequency,  multi-aspect  backscattered  data 
[1,2].  We  have  extended  this  concept  previously  to  the  far-field  antenna 
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Figure  1.  Comparison  of  the  ASAR  and  the  ACSAR  scenario. 


radiation  problem  by  introducing  the  Antenna  SAR  (ASAR)  imaging 
algorithm  [3].  An  ASAR  image  of  a  platform  maps  the  strong  sec¬ 
ondary  radiators  on  the  platform  from  multi-frequency,  multi-aspect 
radiation  data  in  the  far  field.  Such  information  can  be  useful  to  the 
antenna  designer  in  mitigating  platform  effects.  In  this  paper,  we  ex¬ 
tend  the  ASAR  concept  to  the  near-field  antenna  coupling  problem  by 
generating  the  Antenna  Coupling  SAR  (ACSAR)  image  of  the  plat¬ 
form  (see  Fig.  1).  By  collecting  the  multi-frequency,  multi-spatial  cou¬ 
pling  data,  it  is  shown  that  an  ACSAR  image  of  the  platform  can  be 
formed  to  display  the  dominant  scattering  locations  on  the  platform. 
To  achieve  the  required  spatial  diversity,  the  data  are  collected  on  a 

2- D  grid  at  the  receiver  site.  ‘ 

This  paper  is  organized  as  follows.  First,  we  derive  the  ACSAR 
imaging  algorithm.  In  Sec.  2,  it  is  shown  that  under  the  single-bounce 
and  small-bandwidth  approximations,  a  Fourier  transform  relationship 
exists  between  the  multi-frequency,  multi-spatial  coupling  data  and  the 

3- D  locations  of  antenna-platform  interaction.  Hence,  by  3-D  inverse- 
Fourier  transforming  the  coupling  data,  it  is  possible  to  image  these 
scattering  locations  on  the  platform.  This  concept  will  be  demon¬ 
strated  by  using  the  computed  frequency-spatial  data  from  the  Shoot¬ 
ing  and  Bouncing  Ray  (SBR)  technique  [4-6].  In  section  3.  we  present  a 
fast  ACSAR  imaging  algorithm  that  is  specifically  tailored  to  the  SBR 
technique.  By  taking  advantage  of  the  ray  tracing  information  within 
the  SBR  engine,  we  demonstrate  that  a  direct  image-domain  formation 
is  possible  without  resorting  to  multi-frequency,  multi-spatial  calcula¬ 
tions.  The  image-domain  formation  algorithm  can  be  further  acceler- 
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Figure  2.  The  geometry  for  ACSAR  imaging. 


ated  by  the  application  of  the  fast  fourier  transform  (FFT)  algorithm. 
The  fast  algorithm  is  compared  to  the  general  frequency- aperture  al¬ 
gorithm  in  several  numerical  examples.  It  is  also  shown  that  the  total 
computation  time  can  be  reduced  from  hundreds  of  hours  to  minutes 
without  significant  loss  of  fidelity.  In  Sec.  4,  we  further  present  a  sparse 
model  to  represent  ACSAR  imagery.  By  parameterizing  the  ACSAR 
imagery  with  a  set  of  point  radiators  on  the  platform,  it  is  possible  to 
reconstruct  the  3-D  ACSAR  image  and  the  3-D  frequency-spatial  data 
with  a  very  sparse  set  of  radiation  centers. 


2.  ACSAR  IMAGING  ALGORITHM 
2.1  Formulation 

We  shall  first  derive  the  general  ACSAR  imaging  algorithm  that 
utilizes  multi-frequency,  multi-spatial  radiation  data  in  the  near-field 
region  of  the  platform.  We  assume  the  transmitter  and  receiver  setup 
as  depicted  in  Fig.  2.  At  the  receiver  site,  spatial  diversity  on  a  two- 
dimensional  aperture  centered  at  (xo,  0, 0)  is  used  to  achieve  resolution 
in  the  two  cross  range  dimensions.  Similarly,  to  achieve  down  range 
resolution,  frequency  diversity  is  used.  In  addition  to  the  direct  ra¬ 
diation  from  the  transmitter  to  the  receiver,  there  are  contributions 
from  the  antenna-platform  interactions.  The  scattered  electric  field  at 
the  receiver  site  due  to  the  scattering  from  a  point  P(xu  yi,Zi)  on  the 
platform  can  be  written  as 
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Es(k)  =  Ai  ■  e~jkRu  ■  e~jkR2  (1) 

where  A{  is  the  strength  of  the  scattered  field.  Ru  is  the  path  length 
from  the  transmitter  antenna  to  P,  R2  is  the  path  length  from  P 
to  the  receiver  and  k  is  the  free-space  wave  number.  Next,  we  will 
make  two  approximations  to  the  above  equation  to  arrive  at  a  Fourier- 
based  imaging  algorithm.  The  first  assumption,  commonly  used  in 
ISAR  imaging,  is  that  the  radiation  data  are  collected  within  a  certain 
frequency  bandwidth  that  is  small  compared  to  the  center  frequency 
of  operation.  We  will  further  assume  that  the  size  of  the  aperture  at 
the  receiver  site  is  small  compared  to  path  length  R2  .  Combining 
these  assumptions,  we  can  approximate  the  phase  lag  in  the  second 
exponential  term  as  follows: 


kRo  =  kR2i  +  ko(y  cos  ax  +  z  sin  a*  sin  (3X)  (2) 

where  az  and  Pi  are  defined  in  Fig.  2.  is  the  distance  from  P 
to  the  center  of  the  2-D  aperture,  and  the  variables  y  and  z  refer  to 
spatial  locations  in  the  receiver  aperture.  The  scattered  electric  field 
can  thus  be  approximated  by 

Es(k,y,z)  =  Ai  •  e-jklRu+R2i)  ■  €-jk°-ycosai  •  e-J'fco'zsina<sin/?t  (3) 

In  the  above  formula,  a  Fourier  transform  relationship  exists  between 
the  variables  (/c,y,  z)  and  {Rx  =  Ru  4-  R2i,  ux  =  ko  -  cos  a*,  vx  = 
ko  •  sinaisin/?i)  .  By  taking  the  3-D  inverse  Fourier  transform  of  the 
scattered  field  with  respect  to  k ,  y,  and  z,  we  obtain  a  3-D  ACSAR 
image  of  the  platform  as  follows: 

ACSAR(P.  u,  v)  =  IFT  { Es{k ,  y,  z)} 

=  IFT  {  Ai  ■  e~jk  R  ■  e~jko'u  ■  e~jko  V}  (4) 
=  Ai  •  S(R  —  Rx )  •  5(u  —  Ui)  *  S(v  —  Vi) 

Therefore,  by  inverse  Fourier  transforming  the  multi-frequency,  multi- 
spatial  coupling  data,  the  point  scatterer  P  will  manifest  itself  as  a 
peak  in  the  image  at  (RijU^Vi)  with  amplitude  Ai .  In  practice,  since 
the  frequency  bandwidth  and  the  aperture  size  are  not  infinite,  the 
actual  spot  size  of  the  point  scatterer  will  be  inversely  proportional  to 
the  frequency  bandwidth  and  the  aperture  size.  Note  that  we  have  not 
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Figure  3.  The  geometry  for  a  multiple  bounce  mechanism. 

yet  constructed  the  ACSAR  image  in  the  original  (x,  y,  z )  coordinates. 
However,  it  is  straightforward  to  transform  the  ACSAR  image  from  the 
(R,u,v)  to  the  (x,y,  z)  coordinates  by  the  following  transformation 
formulas: 


Xq  -{-  R2  —  2  •  R  *  xo  ■  \/l  +  c 

(5a) 

2  •  (XQ  —  R  ■  y/1  +  c) 

(zo  -  x) 

y  = - 

tan  a  • tan / 3 

(56) 

z  —  tan  / 3  •  (xo  —  x) 

(5c) 

2  1  +  tan2  j 8 

c  —  tan  p  H - ^ - 

tarn  a 

(5  d) 

Since  the  transformation  is  not  linear,  it  is  not  easy  to  estimate  the 
distortion  effect  of  the  transformation  on  the  ACSAR  image.  This 
issue  will  be  addressed  in  the  numerical  examples  as  well  in  Sec  4.3. 
To  summarize,  the  ACSAR  imaging  algorithm  consists  of  three  basic 
steps: 

(1)  Collect  multi-frequency,  multi-spatial  coupling  data  Es(k,y,z) , 

(2)  Take  the  3-D  inverse  Fourier  transform  of  E$(k,  y,  z)  to  form 
ACSAR(jR,tz,  v) , 

(3)  Use  (5)  to  generate  ACSAR (x,y,  z),  the  image  in  the  platform 
coordinates. 

In  the  resulting  ACSAR  image,  the  dominant  scattering  points  on  the 
platform  due  to  platform  interactions  between  the  transmitter  and  the 
receiver  will  be  manifested  as  peaks. 
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Several  comments  are  in  order.  First,  the  present  Fourier-based 
imaging  algorithm  is  derived  based  on  small  bandwidth  and  small 
aperture  approximations.  Since  these  assumptions  are  similar  to  those 
made  in  the  ISAR  imaging  algorithm,  no  attempts  have  been  made  by 
us  in  this  paper  to  quantify  the  errors  related  to  these  approximations. 
Second,  the  present  ACSAR  algorithm  is  based  on  the  single- bounce  as¬ 
sumption.  Therefore,  multiple-bounce  mechanisms  will  not  be  mapped 
correctly  in  the  resulting  ACSAR  image.  This  situation  is  similar  to 
the  standard  ISAR  algorithm.  We  shall  discuss  how  the  multi-bounce 
mechanisms  are  manifested  in  the  ACSAR  image  to  guide  us  in  im¬ 
age  understanding.  We  consider  an  n  -bounce  scattering  mechanism 
as  shown  in  Fig.  3.  The  corresponding  scattered  electric  field  for  this 
mechanism  is  given  by 

Es(k.  y.z)  =  A  •  e”‘-?/c^triPtota,+jR2n^  •  e~Jko-ycosan  .  e-jkO'Zs\nansin0n  ^ 

where  triptotal  =  trip2  +  . . .  +  tripn_1  corresponds  to  the  path  that  the 
wave  travels  from  the  transmitter  to  the  last  hit  point  and  i?2n  is  the 
path  length  from  the  last  hit  point  to  the  center  of  the  aperture  at  the 
receiver  site.  By  taking  the  3-D  inverse-Fourier  transform  of  the  above 
equation,  we  obtain  a  peak  in  the  ACSAR  image  at  ( Ru ,  un.  vn)  where 
Rn  =  triptotal  +  i?2n,  =  ko  •  cos  ocn  and  vn  =  ko  *  sinan  *  sin/?n  . 

We  observe  that  this  corresponds  nearly  to  the  last  hit  point  on  the 
platform.  The  u  and  v  values  (and  the  corresponding  a  and  j3  an¬ 
gles)  correspond  exactly  to  those  of  the  last  hit  points,  while  the  R 
value  corresponds  to  the  cumulative  delay  of  the  mechanism.  This 
means  that  a  multiple-bounce  mechanism  will  be  mapped  in  the  AC¬ 
SAR  image  as  a  point  shifted  from  the  last  hit  point  on  the  platform 
in  a  direction  opposite  to  #2n  ,  the  vector  from  the  last  hit  point  to 
the  receiver.  This  effect  will  be  seen  in  the  example  next.  Finally,  in 
the  ACSAR  algorithm  derivation  only  the  scattered  field  is  considered. 
If  the  total  field  is  applied  to  the  algorithm,  the  direct  radiation  from 
the  transmitter  to  the  receiver  will  be  imaged  at  the  transmitter  loca¬ 
tion.  Since  this  peak  will  likely  be  the  strongest  one  in  the  image,  it 
may  pose  a  dynamic  range  problem  as  the  other  scattering  mechanisms 
due  to  antenna- plat  form  interactions  may  be  overshadowed.  For  sim¬ 
ulation  data,  the  isolation  of  the  scattered  field  contribution  from  the 
direct  antenna  radiation  is  straightforward.  For  measurement  data,  a 
separate  calibration  run  to  measure  the  direct  antenna  radiation  may 
be  required  to  remove  its  contribution  from  the  total  field. 
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Figure  4.  The  geometry  of  the  plate  test  example. 
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Figure  5.  2-D  projected  ACSAR  images  for  the  plate,  (a)  In  (R,  u,  v) 
coordinates,  (b)  In  (x,y,z)  coordinates. 


292 


Ozdemir  and  Ling 


2.2  Numerical  Example 

The  algorithm  presented  above  is  first  tested  using  a  simple  config¬ 
uration  shown  in  Fig.  4.  The  platform  is  a  perfectly  conducting  square 
plate  of  size  10m  by  10m.  The  transmitter  antenna  is  a  half-wave 
dipole  (at  10  GHz)  and  is  located  20m  above  the  plate.  A  total  of  32 
frequencies  are  computed  within  a  0.35  GHz  bandwidth  at  the  center 
frequency  of  10  GHz.  The  simulation  is  carried  out  by  a  modified  ver¬ 
sion  of  the  SBR  code  Apatch  [6].  The  scattered  field  is  collected  on  an 
aperture  centered  at  (20m,  0,  0).  The  collection  aperture  consists  of 
16  x  32  =  512  points,  which  range  from  -0.66m  to  0.56m  in  16  steps 
along  the  y  direction  and  from  -0.46m  to  0.44m  in  32  steps  along  the 
2  direction.  In  the  computed  data  only  the  scattered  field  from  the 
platform  is  considered  and  the  primary  radiation  due  to  the  antenna  is 
not  included.  Using  the  ACSAR  algorithm,  we  first  generate  the  3-D 
ACSAR  images  in  (R,u,v)  domain.  While  applying  the  algorithm,  a 
3-D  Hanning  window  is  used  prior  to  the  FFT  operation  to  suppress 
the  sidelobes  in  the  image.  For  display  purpose,  the  3-D  ACSAR  im¬ 
age  is  projected  onto  the  2-D  R-u ,  R-v  and  u-v  plane  as  shown  in 
Fig.  5(a).  To  view  the  ACSAR  image  in  the  platform  (x,y,  z)  coordi¬ 
nates.  we  perform  the  necessary  transformations  in  (5).  The  resulting 
3-D  ACSAR  image  is  then  projected  onto  the  2-D  y-z,  x-y,  and  x-z 
planes  as  shown  in  Fig.  5(b).  .The  dynamic  range  of  the  images  is  cho¬ 
sen  to  be  35  dB.  According  to  the  geometry,  we  expect  a  peak  in  the 
(x,  y,  z)  image  at  (10m,  -20m ,  0)  corresponding  to  the  specular  point 
on  the  plate.  The  corresponding  R.  u  and  v  values  for  this  point  are 
(44.72m,  186.4,  0).  We  observe  from  Figs.  5(a)  and  5(b)  that  the  peak 
occurs  at  the  expected  location.  However,  we  notice  that  the  specular 
peak  is  not  as  well  focused  in  the  xy  2-plane  as  the  Run-plane .  This 
degradation  is  due  to  the  additional  (R,  u,  u)-to-(x,  y,  2)  transforma¬ 
tion. 

As  the  second  example,  we  used  a  much  more  complex  test  platform 
called  ‘Slicy\  whose  CAD  geometry  is  shown  in  Fig.  6.  Slicy  contains 
a  number  shapes  on  top  of  the  platform  including  a  closed  cylinder, 
an  open  cylinder,  a  corner  reflector  and  a  step  region.  The  platform 
is  assumed  to  be  perfectly  conducting.  A  half-wave  dipole  is  used  as 
the  transmitter  and  is  placed  at  the  origin.  A  total  of  32  frequencies 
are  computed  within  a  0.19  GHz  bandwidth  at  the  center  frequency 
of  10  GHz.  The  scattered  field  is  collected  on  an  aperture  centered 
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Figure  6.  The  geometry  of  the  test  example  Slicy. 


at  (35m,  0,0)  .  There  are  32  x  32  =  1024  spatial  points  on  the  aper¬ 
ture,  ranging  from  -1.14  m  to  1.14  m  along  the  y  direction  and  from 
—0.49  m  to  0.49  m  along  the  z  direction.  The  total  computation  time 
for  the  SBR  simulation  is  about  160  hours  on  an  IBM  RS6000/590 
workstation.  Using  the  ACSAR  algorithm,  we  generate  the  3-D  AC- 
SAR  image  for  the  configuration.  Again,  a  3-D  Hanning  window  is  used 
to  suppress  sidelobes.  The  3-D  ACSAR  image  projected  onto  the  2-D 
R-u,  R-v  and  u-v  planes  are  shown  in  Fig.  7(a).  Fig.  7(b)  shows  the 
projected  ACSAR  image  in  the  platform  (x,  y.  z)  coordinates  along  the 
y-z  ,  x-y  and  x-z  planes.  Also  plotted  in  the  images  is  the  platform 
overlay.  The  dynamic  range  of  the  images  is  35dB.  We  observe  that 
the  dominant  scattering  is  from  the  middle  of  the  platform  around  the 
point  (17.5m,  —40m,  0).  This  is  the  expected  specular  point  on  the 
platform  (labeled  as  mechanism  (i)).  Also  apparent  is  the  scattering 
off  the  curved  region  between  two  steps  (mechanism  (ii)).  This  mecha¬ 
nism  appears  to  be  quite  diffused  in  the  image.  We  believe  this  is  due 
to  the  distortion  effect  from  the  transformation  in  (5).  A  technique 
to  alleviate  this  effect  will  be  discussed  in  Sec.  4.3.  We  also  observe 
some  scattering  mechanisms  from  the  open  cylinder  (mechanism  (iii)). 
Since  the  open  cylinder  is  expected  to  create  multiple  bounces,  the 
contribution  in  the  image  is  shifted  in  the  direction  away  from  the  re- 
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(iii)  u; 


Figure  7.  2-D  projected  ACSAR  images  for  Slicy  (a)  In  ( R,u,v ) 
coordinates,  (b)  In  ( x.y.z )  coordinates. 


ceiver  as  discussed  in  the  last  section.  As  a  result,  it  is  located  around 
the  outside  region  of  the  open  cylinder.  Finally,  the  scattering  off  the 
top  of  the  closed  cylinder  is  labeled  as  mechanism  (iv),  although  it  is 
rather  weak. 

3.  FAST  ACSAR  IMAGE  GENERATION  USING  SBR 

The  ACSAR  image  formation  algorithm  discussed  in  the  last  section 
can  be  applied  to  either  measurement  or  simulation  data.  In  this  sec¬ 
tion,  we  shall  derive  a  fast  ACSAR  imaging  algorithm  specially  tailored 
to  the  SBR  technique.  This  is  accomplished  by  utilizing  the  ray  in¬ 
formation  within  the  SBR  ray-tracing  engine  to  generate  the  ACSAR 
image  without  resorting  to  any  multi-frequency,  multi-spatial  calcula¬ 
tions.  We  first  derive  an  image-domain  formula  of  ACSAR  imaging  [7]. 
Next,  we  implement  a  fast  ray  summation  scheme  [8,9]  based  on  the 
image-domain  formulation  to  reduce  the  total  computation  time.  It  is 
demonstrated  that  an  ACSAR  image  can  be  simulated  at  a  fraction  of 
the  time  as  the  brute-force  frequency-aperture  approach. 
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3.1  Image-Domain  Formulation  of  ACSAR  Imaging 

In  applying  the  SBR  technique  to  the  antenna-platform  interaction 
problem,  rays  are  first  shot  from  the  phase  center  of  the  transmitter 
and  traced  according  to  geometrical  optics.  At  the  exit  point  of  each 
ray  before  they  leave  the  platform  altogether,  a  ray-tube  integration 
is  carried  out  to  find  the  contribution  of  each  ray  to  the  field  at  the 
receiver.  Under  this  construct,  the  scattered  near  field  can  be  written 
as: 

Es(k,  <M)  =  £  ~  Q  •  (&A)exitS(<t>, 0)  ■  e_JkrA  (7) 

i  rays 

where  rA  is  the  position  vector  from  the  last  hit  point  A  to  the 
receiver,  (AA)exit  is  the  cross  section  of  the  exit  ray  tube,  S(<f>,9) 
is  the  normalized  radiation  pattern  from  the  ray  tube  and  C*  is  the 
field  at  the  exit  ray  tube.  For  sufficiently  small  ray  tubes,  5(<£,  6)  can 
be  assumed  to  be  unity.  In  the  image-domain  formulation,  we  will 
assume  that  the  frequency  bandwidth  is  small  compared  to  the  center 
frequency  and  the  aperture  size  at  the  receiver  site  is  small  compared 
to  the  path  length  from  the  last  hit  point  to  the  receiver.  With  these 
assumptions,  we  have  k  •  r  A  =  k  •  R2  +  &o  cos  a  •  y  +  ko  *  sin  a  sin  /?  •  z  . 
We  further  assume  that  the  platform  is  perfectly  conducting.  (For  non- 
perfectly  conducting  platforms,  this  derivation  still  applies  as  long  as 
the  frequency  bandwidth  is  small.)  Using  these  assumptions,  we  can 
let 

(AA)exit  «  <7i  •  e-^tnp)*  (8) 

where  <7*  is  assumed  to  be  independent  of  frequency  and  (trip)*  is  the 
total  path  that  the  i!  th  ray  traveled  from  the  transmitter  to  the  last 
hit  point.  Consequently,  the  scattered  electric  field  at  the  receiver  can 
be  written  as 

Es(k,  y,  z)  =  ^2  *  e~jk{tnp'+R2i)  •  e-ito-eosai-y  .  c-jfco-sinof  sinft-z^9aj 

i  rays 

or  by  the  change  of  variables  similar  to  that  in  Sec.  2,  we  have  simply 
Es(k,  y,  z)  =  J2ai'  e~jkR'  '  e~jUi  V  ‘  e~jVi'Z  (96) 

i  rays 

Substituting  the  above  equation  to  our  ACSAR  imaging  formula,  we 
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get 

ACSAR(ii.  u,  v)  =  lFT3{Es(k,  y,z)} 

=  Y,  Ox  ■  IFT{e^fcK‘}  •  IFT{e-jy'u’}  •  lFT{e~jz  v'} 

i  ravs 

(1°) 

The  inverse-Fourier  transform  operations  in  (10)  can  be  performed  in 
closed  form  with  the  result 

ACSAR(i?.  u.  v)=J2ai-  i2Ak  ■  ejko(R~Ri)  sinc(A k(R  -  J%))} 

i  ravs 

•  {2/coAy  •  sinc(Ay(u  —  u,))} 

*  {2/coAz  •  sinc(Az(t‘  —  v*))} 

(U) 

where  A k  is  the  half  bandwidth  in  frequency.  Ay  and  A z  are  the 
half-lengths  of  the  aperture  in  y  and  z  .  respectively.  We  observe  that 
the  image-domain  formulation  gives  the  contribution  of  each  ray  to  the 
total  ACSAR  image  explicitly.  Since  ,  ( Xi ,  yuzx)  are  available  from 
the  ray  tracing,  it  is  straightforward  to  calculate  the  corresponding 
(Ri,Ui,Vi)  values.  Therefore,  given  (ai.Xi.yi,  Zi)  for  each  ray,  we  can 
form  the  ACSAR  image  directly  by  summing  up  the  weighted  3-D  sine 
functions  in  (11). 

3.2  Fast  Algorithm 

The  image-domain  expression  of  (11)  is  time  consuming  to  carry 
out.  We  next  apply  a  fast  ray  summation  scheme  to  speed  up  the 
calculation  in  (11). We  rewrite  the  image  domain  formula  as 

ACSAR(R,  u,  v)  =  ^  pi  •  h(R  —  R{,  u  -  U{,v  —  Vi)  (12 a) 

i  rays 

where  pi  =  8  •  <7;  •  k% Ak  ■  Ay  ■  A z  and  h(R,  u,  v)  is  the  3-D  ray-spread 
function  given  by 

h(R,  u,v)  =  ejfcofisinc(A k  ■  R)  ■  sinc(Ay  •  u)  ■  sinc(Az  •  v)  (126) 
We  can  cast  (12a)  as  a  convolution  between  a  3-D  weighted  impulse 
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Pi  ■  5{R-  Ri . 

i  rays 


U  ~  Ui ,  V 


*  h(R,  u.v) 


To  speed  up  the  computation  time  of  the  convolution,  we  use  an  FFT- 
based  fast  scheme  proposed  by  Sullivan  [9].  Since  the  3-D  weighted 
impulse  train  g  is  not  uniformly  spaced,  we  first  interpolate  the  func¬ 
tion  onto  a  uniformly  sampled  grid.  To  ensure  the  accuracy  of  com¬ 
putation,  we  oversample  the  grid  at  N  times  the  Nyquist  rate.  For 
the  interpolation,  we  use  a  linear  approximation  to  update  the  closest 
eight  neighbors  depending  on  their  distances  from  the  location  of  the 
original  impulse.  The  detail  of  this  scheme  can  be  found  in  [3]  and 
will  not  be  repeated  here.  After  applying  the  interpolation  scheme, 
the  ACSAR  image  can  be  generated  by 


ACSAR(R,  u,  v)  =  IFFT3(FFT3{gs(R.u,v))-FFT3(h{R,u,v))}  (14) 

where  FFT3  and  IFFT3  are  the  three-dimensional  forward  and  in¬ 
verse  fast  Fourier  transforms,  respectively.  The  final  ACSAR(x,  y,  z) 
image  can  be  generated  using  the  transformations  in  (5). 


3.3  Numerical  Example 


The  fast  approach  presented  above  is  applied  to  the  same  Slicy  ge¬ 
ometry  that  was  used  in  the  frequency-aperture  approach.  In  the  in¬ 
terpolation  grid,  an  over  sampling  ratio  of  N  =  4  is  used  in  all  three 
dimensions  to  ensure  the  accuracy  of  the  computation.  The  results  are 
depicted  in  Figs.  8(a)  and  8(b)  as  2-D  projected  ACSAR  images  in  the 
three  principal  planes.  Fig.  8(a)  shows  the  images  in  the  ( R,u,v )  coor¬ 
dinates,  while  Fig.  8(b)  shows  the  images  in  the  ( x ,  y,  z)  coordinates. 
Again,  a  3-D  Hanning  window  is  applied  and  the  dynamic  range  of 
the  images  is  chosen  to  be  35dB.  By  comparing  these  images  with  the 
images  that  are  generated  by  using  the  frequency-aperture  approach 
(Figs.  7(a)  and  7(b)),  we  observe  a  fairly  good  agreement  where  all 
the  key  features  are  correctly  simulated.  The  rms  error  between  the 
frequency-aspect  result  and  the  image-domain  result  in  the  (R,u,  v) 
domain  is  calculated  to  be  2.13%.  The  image-domain  results  lead  to  a 
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Figure  8.  2-D  projected  ACSAR  images  for  Slicv  using  the  fast  algo¬ 
rithm.  (a)  Image  formed  in  (R.u,v)  coordinates,  (b)  Image  trans¬ 
formed  to  (x,y,z)  coordinates,  (c)  Image  formed  directly  in  (x,j/,z) 
coordinates. 


more  focused  image  due  to  its  small  bandwidth  approximation.  This  is 
consistent  with  our  past  experience  from  ISAR  and  ASAR  algorithms. 
To  show  the  time  savings  that  we  gain  by  using  the  image-domain 
approach,  the  total  computation  times  of  the  frequency-aperture  algo¬ 
rithm  and  the  fast  algorithm  are  compared  in  Table  1.  The  timings 
are  divided  into  the  ray-tracing  time  and  the  ray-summation  time.  As 
it  can  be  seen  from  the  table,  the  ray-tracing  time  of  9  minutes  is  the 
same  for  both  approaches.  Ray  summation  takes  over  160  hours  for 
the  frequency-aperture  approach  and  is  the  dominant  portion  of  the 
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Ray  Tracing 

Ray  Summation 

Total 

Time 

Time 

Time 

Frequency- 

Aperture 

Approach 

9  minutes 

161  hours 

161  hours 

Fast  Algorithm 

9  minutes 

6  minutes 

15  minutes 

Table  1.  ACSAR  simulation  time  for  Slicv  on  an  IBM  RS6000/590 
workstation. 


total  simulation  time.  The  fast  approach  reduces  the  ray  summation 
time  to  only  6  minutes. 

In  addition  to  large  savings  in  simulation  time,  it  is  also  possible  to 
circumvent  the  image  distortion  effect  due  to  the  coordinate  transfor¬ 
mation  by  using  the  image-domain  formulation.  This  is  accomplished 
by  generating  the  ACSAR  image  directly  in  the  platform  (x,  y,  z)  co¬ 
ordinates  in  (12)  instead  of  the  (R,u.v)  coordinates.  To  find  the 
correct  (x^y^Zi)  position  for  a  multi-bounce  ray,  we  move  along  the 
direction  opposite  to  the  vector  from  the  last  hit  point  to  the  receiver. 
The  amount  of  the  shift  is  given  by  the  cumulative  trip  of  the  ray. 
This  has  been  discussed  earlier  in  the  image  interpretation  context 
in  Sec.  2.1.  In  this  manner,  the  image  in  the  platform  coordinate  is 
updated  one  ray  at  a  time  using  the  3-D  ray-spread  function  in  the 
(x,  y,  z)  coordinate  system,  thereby  eliminating  the  need  to  perform 
the  transformation  in  (5).  Fig.  8(c)  shows  the  result  of  generating  the 
image  directly  in  the  (x,y,z)  coordinates.  We  observe  that  all  the 
features  are  now  much  more  focused  than  those  in  Fig.  8(b). 

To  summarize,  we  have  demonstrated  the  feasibility  to  simulate  AC¬ 
SAR  imagery  using  the  SBR  technique  very  rapidly.  This  is  achieved 
by  a  combination  of  the  image-domain  formulation  and  an  FFT-based 
fast  algorithm.  Furthermore,  we  have  shown  that  it  is  possible  to  form 
the  image  in  the  platform  coordinates  directly  by  using  this  approach. 
This  eliminates  the  need  to  perform  the  required  coordinate  transfor¬ 
mation  which  can  defocus  the  features  in  the  image. 


300 


Ozdemir  and  Ling 


4.  A  SPARSE  REPRESENTATION  OF  ACSAR  IMAGERY 
4.1  Point  Radiator  Model  of  ACSAR  Imagery 

In  simulating  the  ACSAR  image  of  the  platform  using  SBR,  tens  of 
thousands  of  rays  are  shot  onto  the  platform.  Yet,  as  we  can  see  from 
the  resulting  data,  the  ACSAR  images  are  in  fact  quite  sparse.  This  is 
because  the  rays  interfere  with  each  other  to  give  rise  to  strong  coherent 
scattering  over  a  small,  localized  area  on  the  platform.  Hence,  it  may 
be  possible  to  accurately  represent  the  ACSAR  image  by  a  finite  set 
of  point  “radiation  centers”  on  the  platform.  This  concept  is  similar 
to  the  scattering  center  representation  in  radar  scattering  problems 
[10,11].  We  shall  adopt  a  point  radiator  model  which  has  the  same 
form  of  (12a)  to  parameterize  the  ACSAR  image 

ACSAR (x,y,z)  =  ^  An  ■  h{x  -  xn,y  -  yn,z  -  zn)  (15) 

n  radiation 
centers 


where  h(x,y,z)  is  given  in  (12b)  and  An  is  the  strength  of  the  point 
radiator  located  at  (xn,yn,zn) .  The  remaining  problem  is  to  deter¬ 
mine  (An,  xn,  yn,  Zn)  to  find  a  radiation  representation  for  the  image. 

4.2  Extraction  Using  CLEAN 


To  perform  the  extraction  process  in  (15),  we  apply  the  image- 
processing  algorithm  CLEAN.  CLEAN  is  a  well-known  extraction  algo¬ 
rithm  first  developed  in  radio  astronomy  [12],  and  has  been  successfully 
applied  for  scattering  center  extraction  [10].  It  is  a  robust  technique 
that  successively  picks  out  the  highest  point  in  the  image  and  removes 
its  point-spread  function  (assumed  to  be  the  ray-spread  function  in  our 
case)  from  the  image.  At  the  n  th  iteration  of  CLEAN,  if  the  strength 
of  the  point  scatterer  is  An  ,  the  3-D  residual  image  can  be  written  as 


[3  —  D  Residual  Image]  n+i  =[3  —  D  Residual  Image]  n 

-  [An  ■  h(x  -xn,y-  yn ,  z  -  zn)\ 


(16) 


The  above  iteration  process  continues  until  the  highest  point  in  the 
residual  image  falls  below  to  a  user-defined  threshold.  Once  we  extract 
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Figure  9.  (a)  The  locations  of  the  150  extracted  radiation  centers,  (b) 
2-D  projected  ACSAR  images  reconstructed  by  using  the  extracted 
radiation  centers. 


the  radiation  centers  from  the  image,  we  can  reconstruct  the  multi¬ 
frequency  and  multi-spatial  data  by  the  formula: 

Es(k,  y,  z)  =  An-  e~ikRn  ■  e~jUny  ■  e~^"'2  (17) 

n  radiation 
centers 


where  (Rn,un,vn)  are  related  to  ( xn,yn,zn )  through  the  transforma¬ 
tion  described  earlier. 

We  apply  the  CLEAN  extraction  algorithm  to  the  ‘Slicy:  example. 
We  extract  a  total  of  150  point  radiators  from  the  3-D  ACSAR  image 
shown  in  Fig.  8(c).  The  locations  of  the  extracted  radiation  centers  are 
projected  on  three  different  2-D  planes  and  are  shown  as  small  circles 
in  Fig.  9(a).  As  it  can  be  seen  from  the  figures,  most  of  the  radiation 
centers  are  located  on  the  platform  at  y  =  —40m .  The  remainders  are 
located  on  the  top  of  tall  cylinder  and  around  the  step  region.  We  then 
reconstruct  the  ACSAR  image  using  the  radiation  centers  and  they  are 
shown  in  Fig.  9(b).  We  use  the  same  dynamic  range  of  35dB  when 
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- Reconstructed 


Figure  10.  Comparison  of  the  original  and  the  reconstructed  field 
patterns,  (a)  Frequency  sweep  at  (35,  0,  0),  (b)  Spatial  comparison 
at  10  GHz  along  (35.  y ,  0),  (c)  Spatial  comparison  at  10  GHz  along 
(35,  0,  z). 


plotting  the  images.  By  comparing  the  reconstructed  ACSAR  images 
to  the  reference  image  in  Fig.  8(c),  we  see  that  the  reconstruction 
is  good.  The  rms  error  between  the  original  and  the  reconstructed 
ACSAR  images  is  calculated  to  be  less  than  0.1%.  Next,  we  reconstruct 
the  frequency  and  the  spatial  data  by  using  the  formula  in  (17).  The 
original  (solid)  and  the  reconstructed  (dashed)  patterns  are  shown  in 
Fig.  10.  In  Fig.  10(a),  the  frequency  sweep  comparison  is  made  at 
the  center  point  of  the  receiver.  In  Fig.  10(b)  and  10(c),  the  spatial 
variation  comparisons  are  made  at  the  10  GHz  center  frequency  along 
the  y-axis  and  along  the  z-axis ,  respectively.  We  observe  a  fairly 
good  agreement  for  the  frequency-spatial  reconstruction  by  using  the 
radiation  center  model.  Therefore,  we  can  represent  complex  coupling 


ACSAR-antenna  coupling  synthetic  aperture  radar  imaging 


303 


-10 


10 


0  10  _2j  30 


0 


10  30 


Figure  11.  2-D  projected  ACSAR  images  based  on  frequency-aperture 
data,  formed  by  first  parameterizing  the  ( R ,  u ,  v)  image  using  CLEAN 
and  reconstructing  the  image  in  the  (x,y,z)  platform  coordinates. 


data  between  antennas  on  a  platform  by  using  only  a  small  set  of 
point  radiators.  Such  sparse  data  representation  may  be  useful  for 
data  compression  and  system  simulation  applications. 

4.3.  Eliminating  Coordinate  Transform  Distortion  Using 
CLEAN 

We  now  describe  a  method  that  takes  advantage  of  the  CLEAN  al¬ 
gorithm  to  eliminate  the  image  distortion  resulting  from  the  ( R ,  u,  v)  - 
to-  (x,  y7  z)  coordinate  transformation  in  the  frequency-aperture  image 
formation  algorithm  described  in  Sec.  2.  In  Sec.  3.3,  we  have  already 
shown  that  this  can  be  accomplished  by  using  the  image-domain  for¬ 
mulation.  However,  the  approach  is  only  applicable  when  coupled  to 
the  SBR  simulation.  When  data  from  measurement  or  other  types  of 
computational  electromagnetics  simulators  are  used,  we  must  adopt  a 
different  strategy.  Our  approach  is  to  first  parameterize  the  image  in 
the  (R,u7  v)  coordinate  system  using  the  CLEAN  algorithm.  Once  the 
image  is  parameterized  by  radiation  centers  located  at  (Rn,  un,  vn) ,  we 
use  (5)  to  find  the  corresponding  (xn,  yn,  zn)  location  of  each  radia¬ 
tion  center  in  the  platform  coordinate.  Then,  the  construction  of  the 
image  in  the  platform  coordinate  is  carried  out  using  (15).  Note  that 
the  main  advantage  of  this  approach  is  that  the  point-spread  response 
in  the  final  platform  image  is  well-controlled.  Consequently,  the  dis¬ 
torted  point-spread  response  in  the  (x,y,  z)  coordinates  is  avoided  in 
the  platform  image. 
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As  an  example,  we  apply  this  methodology  to  the  ACSAR  image  in 
Fig.  7(a)  by  extracting  a  total  of  150  radiation  centers.  After  trans¬ 
forming  these  radiation  centers  into  the  platform  coordinates,  we  form 
the  ACSAR  image  directly  in  the  (x,y,z)  coordinates.  The  result 
is  shown  in  Fig.  11.  By  comparing  the  new  images  to  the  images  in 
Fig.  7(b).  we  observe  that  the  mechanisms  are  much  better  focused 
using  this  approach.  It  is  therefore  the  preferred  step  to  generate  the 
final  platform  image  using  the  general  frequency-aperture  algorithm. 

5.  CONCLUSION 

The  ACSAR  (antenna  coupling  synthetic  aperture  radar)  imaging  al¬ 
gorithm  has  been  introduced  to  map  platform  scattering  locations  in 
antenna  coupling  scenarios.  It  is  shown  that  under  the  single-bounce 
and  small-bandwidth  approximations,  a  Fourier  transform  relationship 
exists  between  the  multi-frequency,  multi-spatial  coupling  data  and  the 
3-D  locations  of  antenna-platform  interaction.  Hence,  by  3-D  inverse- 
Fourier  transforming  the  coupling  data,  it  is  possible  to  image  the  scat¬ 
tering  locations  on  the  platform.  This  concept  has  been  demonstrated 
by  using  the  computed  frequency-spatial  data  from  the  SBR  technique. 
However,  the  algorithm  is  general  and  can  be  applied  to  other  simu¬ 
lation  data  or  even  measurements.  Furthermore,  we  have  presented  a 
fast  ACSAR  imaging  algorithm  that  is  specifically  tailored  to  the  SBR 
technique.  By  taking  advantage  of  the  ray  tracing  information  within 
the  SBR  engine,  we  have  demonstrated  that  a  direct  image-domain  for¬ 
mation  is  possible  without  resorting  to  multi-frequency,  multi-spatial 
calculations.  The  image-domain  formation  algorithm  can  be  further  ac¬ 
celerated  by  the  application  of  the  FFT  algorithm.  The  fast  algorithm 
has  been  shown  to  reduce  the  total  simulation  time  from  hundreds  of 
hours  to  minutes  without  significant  loss  of  fidelity.  Finally,  we  have 
presented  a  method  to  extract  a  sparse  point-radiator  model  to  rep¬ 
resent  ACSAR  imagery.  This  is  accomplished  by  parameterizing  the 
ACSAR  image  with  a  set  of  point  radiators  on  the  platform  using  the 
CLEAN  algorithm.  Once  such  a  sparse  set  of  radiation  centers  have 
been  extracted,  it  is  possible  to  reconstruct  the  3-D  ACSAR  image  and 
the  3-D  frequency-spatial  data  very  rapidly  with  good  fidelity. 
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ABSTRACT 

An  adaptive  approach  is  proposed  to  construct  ISAR  images  from  unevenly 
undersampled  data  in  the  angular  domain.  The  algorithm  uses  an  adaptive  scattering 
feature  extraction  engine  in  place  of  the  Fourier  transform  in  the  image  construction 
procedure.  The  algorithm  entails  searching  and  extracting  out  individual  target  scattering 
features  one  at  a  time  in  an  iterative  fashion.  Therefore  the  strong  interference  between 
different  target  scattering  features  can  be  reduced.  After  all  the  main  features  are 
extracted,  they  can  be  displayed  in  the  same  image  plane  without  the  aliasing  effect  in  the 
final  ISAR  image.  The  algorithm  is  verified  by  constructing  the  ISAR  image  from  the 
simulated  point  scatterer  data  as  well  as  the  chamber  measurement  data  of  the  scaled 


VFY-218  airplane. 


I.  Introduction 


Inverse  Synthetic  Aperture  Radar  (ISAR)  imaging  is  a  well-known  tool  for 
pinpointing  target  scattering  features  for  signature  diagnostic  and  target  identification 
purposes.  Constructing  an  ISAR  image  requires  data  collection  in  both  the  frequency  and 
angular  dimensions.  If  the  data  are  evenly  sampled  and  the  sampling  rate  is  dense  enough, 
an  ISAR  image  can  be  obtained  by  using  a  two-dimensional  Fourier  transform  algorithm 
[1].  In  this  paper,  we  address  the  case  when  the  angular  data  are  unevenly  undersampled. 
This  scenario  may  arise  in  real-world  data  collection  when  the  target  is  fast  maneuvering 
with  respect  to  the  radar  pulse  repetition  interval  so  that  the  angular  look  on  the  target  by 
the  radar  is  not  dense  enough  to  satisfy  the  Nyquist  sampling  rate.  This  will  cause  serious 
aliasing  error  in  the  cross  range  dimension  if  the  regular  Fourier  transform  is  used.  We 
propose  an  algorithm  to  overcome  the  aliasing  effect  in  the  cross  range  dimension  and 
construct  ISAR  images  from  undersampled  data.  This  algorithm  can  also  be  easily 
extended  to  deal  with  undersampling  in  both  the  frequency  and  angular  domains. 

The  algorithm  uses  an  adaptive  scattering  feature  extraction  engine  in  place  of  the 
Fourier  transform  in  the  image  construction  process.  The  original  concept  of  adaptive 
feature  extraction  was  proposed  in  [2,3].  It  has  been  applied  to  ISAR  image  processing  in 
the  joint  time-frequency  space  for  resonant  scattering  mechanism  extraction  [4],  target 
motion  compensation  [5]  and  Doppler  interference  removal  [6].  In  contrast  to  the  Fourier 
transform  where  the  signal  is  projected  to  all  the  image  domain  bases  simultaneously,  the 
adaptive  algorithm  searches  and  extracts  the  individual  target  scattering  features  one  at  a 
time  in  an  iterative  fashion.  When  applied  to  the  present  problem,  the  severe  aliasing 
error  caused  by  the  interference  between  different  target  scattering  features  can  be 


avoided.  Therefore,  after  all  the  main  features  are  extracted,  they  can  be  displayed  in  the 
image  plane  without  the  aliasing  effect  in  the  final  ISAR  image.  This  paper  is  organized 
as  follows.  In  Sec.  II,  we  describe  the  adaptive  feature  extraction  algorithm  based  on  a 
point  scatterer  model.  In  Sec.  HI,  we  verify  the  algorithm  by  using  simulated  data 
consisting  of  point  scatterers.  In  Sec.  IV,  the  algorithm  is  further  tested  using  the  chamber 
measurement  data  of  the  model  VFY-218  airplane  [7].  It  is  found  that  a  reasonable  ISAR 
image  can  be  constructed  from  seriously  undersampled  data.  Finally,  discussions  on  the 
limitation  of  the  algorithm  and  some  conclusions  are  given  in  Sec.  V. 


II.  Adaptive  Feature  Extraction  (AFE)  Algorithm 
In  standard  ISAR  image  construction,  the  target  is  assumed  to  be  a  collection  of 
point  scattering  centers.  Under  the  small-angle  approximation,  the  scattered  field  from 
the  target  observed  versus  frequency  and  angle  can  be  written  as: 

E(f,  0)  =  f^OiXi ,  yi  )e~2jkx'  <*x0-2jkyi  sine  s  £0(JC.  f  y.  )e-'ljkxi-2jkc0yi  (}) 
i=l  i=l 


where  0(Xi  ,yj  is  the  amplitude  of  the  ith  scattering  center,  k  is  the  free  space  wave 
number  and  kc  corresponds  to  the  wave  number  at  the  center  frequency,  x,  and  y,  denote 
the  down  range  and  cross  range  dimensions,  respectively.  We  shall  assume  that  the 
sampling  in  frequency  is  evenly  spaced  and  dense  enough  to  satisfy  the  Nyquist  criterion. 
This  is  a  reasonable  assumption  since  the  frequency  parameter  is  completely  controlled 
by  the  radar.  Thus  the  range  profile  versus  angle  data  can  be  generated  from  the 
frequency-aspect  data  by  applying  a  1-D  Fourier  transform  (which  can  be  implemented 


using  the  FFT  algorithm)  along  the  frequency  dimension.  We  shall  denote  the  result  as 

R(x,d): 


N 


R(x, 6)  =  X 0(Xi ,  yi )Sx{x- xt  )e 


-2jkcyi 


i= 1 


(2) 


In  the  above  expression,  Sx(x-Xi)  is  the  point  spread  function  due  to  the  finite-length 
frequency  domain  data.  For  example,  for  unwindowed  frequency  data  the  point  spread 
function  in  down  range  is  given  by: 


Ak 

Sx  =  Ak  ■  sinc[ — (x  -  )]  •  exp(-./2A:cx,- ) 

2 


(3) 


where  Ak  is  the  frequency  bandwidth.  Similarly,  the  cross  range  information  can  also  be 
obtained  from  angular  data  via  a  1-D  Fourier  transform  of  R(x  ,6)  along  the  angular 
dimension.  The  resulting  image  l(x  ,y)  is: 


I(x  ,y  )  =  j* R(x  ,d)e2jk'y°d0  =  '£O(xi,yi)Sx(x-xi)\e2jk'9(y  yi)d6 

(=1 


N 

=  X  °(Xi '  yi  ">SX  (x  -  xi  )Sy{y-yi) 
1=1 


(4a) 


where 

Sy(y-yi)  =  \e2ik‘S(-y-y<)de  (4b) 

is  the  cross-range  point  spread  function  due  to  the  finite-length  angular  domain  data.  In 
the  ideal  case,  the  target  rotates  at  a  constant  velocity  so  that  the  observation  angle  is 
proportional  to  the  dwell  time  and  is  sampled  evenly.  Thus  the  integral  in  (4b)  can  be 
calculated  using  an  FFT  algorithm.  However,  if  the  data  are  not  evenly  sampled  in  angle, 
the  integration  must  be  carried  out  numerically  in  place  of  the  FFT  algorithm.  We  shall 


now  turn  our  attention  to  the  undersampling  issue.  If  the  data  are  sampled  densely 
enough  so  that  the  numerical  integration  can  be  carried  out  accurately,  the  integration 
result  Sy  should  be  a  well  localized  function  with  its  peak  at  y,  while  rapidly  decaying 
away  from  the  peak,  similar  to  the  point  spread  function  Sx.  Under  this  condition,  the 
resulting  image  I(x,  y )  will  be  a  clean  image  with  good  indication  of  the  amplitudes  and 
positions  of  the  target  point  scattering  features.  However,  when  the  data  are 
undersampled,  the  numerical  integration  in  (4b)  will  result  in  large  aliasing  error  that 
shows  up  as  high  sidelobes  in  Sy.  Consequently,  the  constructed  image  will  contain 
strong  interference  between  the  scattering  features.  This  effect  can  be  interpreted  as  the 
loss  of  orthogonality  of  the  Fourier  bases  under  the  undersampled  condition.  For  the 
scattering  signal  from  a  given  target,  the  Nyquist  sampling  criterion  requires  the 
maximum  interval  between  two  consecutive  angles  be  smaller  than  where  is 

the  maximum  cross  range  dimension  of  the  target.  When  this  criterion  is  met,  the  sidelobe 
noise  is  outside  the  target  region,  and  (4)  gives  the  correct  image  of  the  target.  Otherwise, 
the  sidelobe  noise  will  overlap  with  the  target  features,  thus  causing  aliasing  and 
corrupting  the  image. 

In  the  proposed  approach,  we  use  an  adaptive  feature  extraction  algorithm  in  place 
of  the  Fourier  processing.  To  avoid  the  sidelobe  interference  between  different  image 
bases,  an  iterative  procedure  is  employed.  Instead  of  projecting  the  signal  onto  all  the 
exponential  bases  simultaneously,  we  search  out  the  strongest  point  scattering  feature  in 
the  cross-range  domain  and  remove  it  from  the  original  signal.  Then  the  search  is 
repeated  for  the  remainder  signal  and  the  point  scattering  features  are  extracted  one  at  a 
time  until  the  energy  of  the  residue  signal  is  smaller  than  a  preset  threshold.  The  search 


procedure  is  carried  out  by  calculating  the  integral  in  (4b)  for  all  the  points  in  cross  range 
but  saving  only  the  maximum  value  and  position,  i.e., 

[Bp,yp]  =  max[Ip(x,y)]  (5) 

yP 

where  p  denotes  the  pth  stage  of  the  iterative  procedure.  The  remainder  signal  is 
produced  by  subtracting  out  the  pth  feature: 

Rp+l  (x,  6)  =  Rp  (x,  0)  -  Bpe~2jkcyP6  (6) 

The  convergence  of  the  above  procedure  is  guaranteed  and  the  mathematical  proof  is 
given  in  [2].  The  advantage  of  such  an  iterative  procedure  is  that  each  time  we  extract  out 
the  strongest  feature,  we  also  eliminate  its  interference  on  the  other  features.  It  should  be 
noted  that  unevenly  sampling  is  a  prerequisite  to  ensure  that  there  is  no  ambiguity  in  the 
strongest  features,  since  evenly  undersampled  data  will  result  in  multiple  positions  of  the 
strongest  features.  After  all  the  features  are  extracted  out,  we  can  construct  an  ISAR 
image  using  the  amplitudes  and  positions  of  the  point  scatterers. 

III.  Test  Results  Using  Simulated  Data 

To  verify  the  algorithm,  we  first  simulate  scattering  data  using  1 1  point  scattering 
centers  with  different  amplitudes  and  positions.  The  radar  center  frequency  is  assumed  to 
be  10.4  GHz  and  the  bandwidth  is  1  GHz.  The  total  angular  observation  aperture  is  5°. 
The  frequency  sampling  is  fixed  at  256  points.  The  sampling  in  the  angular  domain  is 
varied  to  compare  the  performance  of  the  algorithm.  Shown  in  Fig.  1(a)  is  the  ISAR  image 
constructed  using  the  Fourier  transform  from  128  evenly  sampled  points  in  the  angular 
domain.  No  aliasing  error  occurs  in  the  image  because  the  angular  sampling  satisfies  the 


Nyquist  rate.  However,  if  the  same  128  points  are  not  evenly  sampled  within  this  fixed 
angular  range,  the  sampling  no  longer  satisfies  the  Nyquist  rate  because  some  of  the 
intervals  between  two  consecutive  angles  is  greater  than  n/kcymax-  Fig.  1(b)  is  the  ISAR 
image  constructed  using  the  Fourier  transform  from  128  randomly  sampled  angular  data 
points.  The  aliasing  effect  can  be  clearly  seen  from  Fig.  1(b).  It  behaves  like  sidelobe 
noise  smearing  across  the  cross  range  dimension.  The  aliasing  error  gets  more  severe 
when  fewer  sampled  points  are  used.  This  is  illustrated  in  Fig.  2(a),  which  is  the  image 
obtained  using  the  Fourier  transform  from  64  randomly  sampled  points  in  the  angular 
domain.  We  next  apply  the  AFE  algorithm  to  the  same  set  of  undersampled  data.  The 
resulting  image  is  plotted  in  Fig. 2(b).  Comparing  Fig.  2(b)  with  the  reference  image  in 
Fig.  1(a),  we  can  see  that  the  aliasing  error  is  removed  and  the  amplitudes  and  positions 
of  all  the  scattering  center  features  are  well  restored  in  the  image  through  AFE.  The 
performance  of  the  AFE  is  further  demonstrated  when  the  angular  sampling  is  lowered  to 
32  random  points.  The  ISAR  image  constructed  using  the  Fourier  transform  and  the  AFE 
algorithm  are  plotted  respectively  in  Figs.  3(a)  and  3(b).  We  observe  that  AFE  is  still 
able  to  remove  the  aliasing  error  in  the  image  and  restores  all  the  point  scattering  features 
while  the  image  obtained  using  the  Fourier  transform  is  badly  overwhelmed  by  aliasing 
error. 


IV.  Test  Results  Using  Measurement  Data 

To  examine  the  applicability  of  the  algorithm  on  real  target  scattering  data,  we 
reconstruct  the  radar  image  of  a  model  (1:30  scale)  VFY-218  airplane  using 
undersampled  chamber  measurement  data  [7].  The  original  measurement  data  consist  of 


201  samples  in  a  40-degree  azimuth  aperture,  and  401  samples  in  frequency  from  8  GHz 
to  16  GHz.  To  construct  an  ISAR  image,  we  first  polar-reformat  the  frequency-aspect 
data  to  the  (Kx,  Ky)  space,  since  the  observation  angular  window  involved  is  too  large  to 
use  the  small-angle  approximation  in  (1).  (However,  for  full-size  targets  the  required 
angular  swath  is  typically  small  enough  so  that  polar  reformatting  is  not  very  important.) 
The  reformatted  data  consist  of  401  samples  in  Kx  and  438  samples  in  Ky.  For  the 
observation  azimuth  aperture  from  110  degrees  to  150  degrees,  the  ISAR  image  is  first 
generated  by  FFT  to  the  X-Y  plane  and  is  shown  with  the  airplane  overlay  in  Fig.  4.  The 
point  scattering  features  can  clearly  be  seen.  Next  we  test  our  algorithm  by  generating  an 
undersampled  data  set  in  Ky.  This  is  approximately  the  same  as  undersampling  in  angle. 
(Again,  for  full  size  targets,  this  approximation  gets  better.)  Note  that  the  original 
measurement  data  set  is  oversampled.  We  can  downsample  Ky  by  at  most  12  times 
without  violating  the  Nyquist  sampling  rate,  which  means  at  least  36  points  must  be 
evenly  sampled  in  Ky.  The  test  is  conducted  as  follows.  First  we  randomly  choose  24  out 
of  the  438  points  (in  this  example,  at  [1,  27,  30,  36,  77,  79,  1 10,  143,  169,  199,  211,  260, 
272,  280,  289,  320,  329,  333,  367,  378,  390,  420,  438]).  Note  that  the  maximum 
sampling  interval  is  49  points  which  is  much  larger  than  the  Nyquist  sampling  of  every  12 
points.  Serious  aliasing  occur  if  the  regular  Fourier  transform  algorithm  is  used,  as  shown 
by  the  ISAR  image  displayed  in  Fig.  5(a).  All  the  features  are  overlapped  with  sidelobe 
noise  so  that  no  point  scattering  features  on  the  airplane  can  be  distinguished  anymore. 
Next,  we  apply  the  AFE  algorithm  for  each  range  cell  of  the  image.  The  image  is 
reconstructed  and  shown  in  Fig.  5(b).  Comparing  Fig.  5(b)  with  Fig.  4,  we  can  see  the 
main  features  of  the  airplane  in  Fig.  4  are  all  well  reconstructed  in  Fig.  5(b).  We  do 


observe  some  noisy  spots  outside  the  target  at  the  lower  dynamic  ranges  in  Fig.  5(b). 
This  low-level  noise  occurs  at  about  25dB  down  from  the  key  features  and  presents  a 
dynamic  range  limitation  of  the  present  AFE  algorithm.  We  also  notice  that  the  weak 
dispersive  cloud  over  the  right  wing  of  the  airplane  due  to  the  exhaust  duct  in  Fig.  4  does 
not  appear  fully  in  Fig.  5(b).  The  possible  reason  is  that  the  algorithm  assumes  a  point- 
scatterer  basis,  so  dispersive  returns  may  not  be  well  reconstructed.  The  same  test  has 
been  done  for  another  set  of  observation  angles.  This  time  the  radar  observation  angle  is 
in  the  front  of  the  airplane  from  10  degrees  to  50  degrees.  The  original  ISAR  image  for 
this  angle  is  plotted  in  Fig.  6.  The  images  reconstructed  from  the  undersampled  data  using 
the  Fourier  processing  and  the  AFE  are  plotted  respectively  in  Figs.  7(a)  and  7(b).  Again 
we  observe  a  good  reconstruction  of  the  key  scattering  features  in  Fig.  7(b). 

V.  Conclusions  and  Discussions 

An  adaptive  feature  extraction  algorithm  has  been  proposed  in  place  of  regular 
Fourier  processing  in  ISAR  image  construction  from  unevenly  undersampled  data.  Based 
on  a  point  scattering  model  of  the  target,  the  Nyquist  sampling  criterion  could  be 
overcome  using  this  adaptive  algorithm.  The  algorithm  has  been  successfully  tested  using 
both  simulated  data  and  chamber  measurement  data.  The  images  reconstructed  from 
unevenly  undersampled  data  agree  very  well  with  the  original  image  without  aliasing 
error.  The  algorithm  has  also  been  tested  on  measured  data  from  in-flight  targets  with 
good  success. 

Three  remarks  are  in  order.  First,  although  we  have  focused  on  the  angular 
undersampling  issue,  the  AFE  algorithm  can  be  easily  extended  to  two  dimensions  to 


process  undersampled  data  in  both  the  frequency  and  angular  domains.  Second,  we  expect 
the  performance  of  the  AFE  algorithm  to  be  strongly  dependent  on  how  well  the  model 
used  reflects  the  real  physical  phenomenon.  In  this  paper,  the  point  scatterer  model  is 
assumed  for  AFE,  which  is  a  coarse  approximation  of  real  target  scattering.  It  is  also 
possible  to  apply  AFE  based  on  more  sophisticated  models  such  as  point  scatterers  with 
frequency  dependence.  In  that  case,  the  basis  functions  in  AFE  can  be  extended  to 
include  a  broader  dictionary.  The  tradeoff  is  expected  to  be  better  performance  at  the 
expense  of  algorithm  complexity.  Third,  as  we  have  observed,  the  dynamic  range  the 
AFE  algorithm  can  achieve  is  limited.  The  reason  is  that  the  manner  in  which  the 
scattering  features  are  selected  in  the  algorithm  is  not  without  error.  Both  the  amplitude 
and  position  of  the  strongest  basis  are  in  fact  contaminated  by  sidelobes  of  the  other 
weaker  bases.  Possible  ways  to  improve  the  dynamic  range,  such  as  incorporating  a 
relaxation  procedure  to  modify  the  previously  found  scattering  features  [8],  are  currently 
being  investigated. 
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Fig.l.  (a)  Simulated  ISAR  image  from  128  points  evenly  sampled  in  the  angular 
domain  using  Fourier  transform,  (b)  Simulated  ISAR  image  from  128  points 
randomly  sampled  in  the  angular  domain  using  Fourier  transform. 
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Fig.2.  (a)  Simulated  ISAR  image  from  64  points  randomly  sampled  in  the  angular 
domain  using  Fourier  transform.  (b)Simulated  ISAR  image  from  64  points 
randomly  sampled  in  the  angular  domain  using  the  AFE  algorithm. 
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Fig. 3.  (a)  Simulated  ISAR  image  from  32  points  randomly  sampled  in  the  angular 
domain  using  Fourier  transform.  (b)Simulated  ISAR  image  from  32  points 
randomly  sampled  in  the  angular  domain  using  the  APE  algorithm. 
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Fig.4  The  ISAR  image  constructed  at  130  degree  azimuth  from  the  original  data 
using  FFT . 
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Fig.5(a)  The  ISAR  image  constructed  at  130  degree  azimuth  from  randomly 
undersampled  data  using  Fourier  transform. 
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Fig.5(b)  The  ISAR  image  constructed  at  130  degree  azimuth  from 
randomly  undersampled  data  using  the  AFE  algorithm. 
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Fig.6  The  ISAR  image  constructed  at  30  degree  azimuth  from  the  original  data 
using  FFT . 
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Fig. 7(a)  The  ISAR  image  constructed  at  30  degree  azimuth  from  randomly 
undersampled  data  using  Fourier  transform. 
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Fig.7(b)  The  ISAR  image  constructed  at  30  degree  azimuth  from  randomly 
undersampled  data  using  the  APE  algorithm. 
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Abstract 

A  methodology  based  on  the  genetic  algorithm  (GA)  is  proposed  to  determine  the 
equivalent  impedance  boundary  condition  (IBC)  for  corrugated  material  coating 
structures.  In  this  approach,  rigorous  solutions  of  the  reflection  coefficients  at  a  number 
of  incident  angles  are  first  calculated  using  a  periodic  method  of  moments  (MoM)  solver. 
The  IBC  model  is  used  to  predict  the  reflection  coefficients  at  the  same  observation 
angles.  The  model  coefficients  are  then  optimized  using  GA  so  that  the  difference 
between  the  approximated  and  the  MoM-predicted  reflection  coefficients  is  minimized. 
The  genetic  algorithm  proves  efficient  in  obtaining  an  optimal  IBC  model.  The  resulting 
IBC  model  can  be  readily  incorporated  into  an  existing  computational  electromagnetics 
code  to  assess  the  performance  of  the  corrugated  coating  when  mounted  on  complex 
platforms. 


I.  Introduction 


It  is  well  known  that  the  impedance  boundary  condition  (EBC)  approximation  is  an 
efficient  way  to  model  complex  structures  such  as  material  coatings  and  sub-skinline 
features  [1-3].  It  replaces  the  original  volumetric  structure  with  a  surface  impedance  so 
that  the  problem  dimension  is  reduced  by  one.  Thus,  large  savings  in  computational 
resources  can  be  achieved  in  the  analysis  of  the  original  problem.  However,  to  determine 
a  simple  EBC  for  an  arbitrary  structure  that  is  valid  over  a  wide  range  of  incident  angles, 
polarizations  and  frequencies  is  a  non-trivial  task.  In  this  paper,  we  set  out  to  develop  a 
methodology  to  determine  the  equivalent  EBC  model  for  a  corrugated  coating  structure 
backed  by  a  conducting  surface  (see  Fig.  1).  The  corrugation  of  the  surface  is  assumed  to 
be  periodic  in  one  dimension  along  the  x-axis.  Of  interest  is  an  IBC  model  that  is  valid 
over  a  large  range  of  incident  angles  in  both  the  6  and  0  directions.  Our  objective  is  to 
establish  a  robust  methodology  such  that  the  resulting  EBC  model  can  be  used  in  place  of 
the  actual  coating  structure  in  subsequent  analysis  and  design  involving  complex 
platforms. 

The  problem  at  hand  is  difficult  since  the  scattering  characteristics  of  the  corrugated 
surface  is  strongly  dependent  on  polarization  and  incident  angle.  The  standard  EBC  used 
for  flat  coatings  accounts  for  neither  the  aniostropic  nor  the  angular  behavior  of  the 
scattering  characteristics  of  the  corrugated  surface.  Some  improved  impedance  boundary 
conditions  have  been  proposed  in  the  literature,  including  the  tensor  impedance  boundary 
condition  (TEBC)  [1]  and  the  generalized  impedance  boundary  condition  (GEBC)  [2], [4]. 
TEBC  usually  works  only  for  a  very  limited  range  of  incident  angles.  GEBC  improves  the 
accuracy  of  the  IBC  model  by  including  higher  order  derivatives  of  the  fields  on  the 


surface.  However,  it  cannot  be  easily  implemented  in  existing  MoM  solvers  since  it 
requires  further  reformulation  in  the  integral  equation.  A  resistive  boundary  condition 
(RBC)  has  been  reported  that  works  well  over  large  incident  angles  for  2-D  planar 
periodic  surfaces  [5].  However,  it  is  limited  to  surfaces  with  very  small  periods. 
Furthermore,  the  choice  for  the  position  of  the  equivalent  impedance  surface  is  not 
obvious  for  the  corrugated  structure. 

Our  proposed  approach  to  this  problem  is  based  on  the  genetic  algorithm  (GA).  First, 
we  compute  the  reflection  coefficients  from  the  corrugated  surface  over  a  number  of 
incident  angles  and  polarizations  using  a  periodic  method  of  moments  (MoM)  solver  [6]. 
The  resulting  reflection  coefficients  constitute  our  reference  data  base.  Next,  a  simple 
periodic  EBC  model  is  proposed,  from  which  we  can  derive  an  expression  for  the 
reflection  coefficients.  In  the  GA  step,  we  optimize  the  model  coefficients  so  that  the 
difference  between  the  IBC-predicted  and  the  MoM-predicted  reflection  coefficients  is 
minimized.  GA  searches  the  entire  parameter  space  in  a  way  similar  to  natural  evolution 
and  arrives,  after  many  generations,  at  the  best  parameters  for  the  model. 

This  paper  is  organized  as  follows.  In  Section  II,  the  MoM  solution  of  the  problem, 
the  IBC  model  formulation  and  the  GA  optimization  are  discussed  as  the  steps  of  the  EBC 
determination.  Numerical  results  are  provided  in  Section  III  to  verify  the  effectiveness  of 
the  approach  in  a  number  of  corrugated  geometries.  A  3-D  scattering  example  is  also 
given  to  demonstrate  the  utility  of  the  resulting  EBC  model. 


II.  Methodology  for  Determining  the  Optimal  IBC  Model 


In  this  section,  we  describe  the  proposed  methodology  for  determining  the  equivalent 
IBC  of  a  corrugated  coating  using  the  genetic  algorithm.  In  the  first  step,  the  reflection 
coefficients  from  the  coating  are  computed  using  the  MoM  at  multiple  incident  angles  to 
serve  as  the  reference  data  of  the  model.  The  MoM  solution  for  the  corrugated  coating 
structure  in  Fig.  1  has  been  formulated  earlier  in  [6].  The  formulation  entails  dividing  one 
cell  of  the  grating  into  different  homogeneous  regions  according  to  the  material  layers  as 
shown  in  Fig.  2  (a).  Boundary  integral  equations  are  first  obtained  for  each  region.  Field 
continuity  at  region  interfaces  and  periodic  boundary  conditions  at  cell  boundaries  are 
then  enforced.  The  fields  in  the  top  half-space  are  expanded  into  a  sum  of  Floquet 
harmonics  and  are  matched  to  the  fields  in  the  lower  region  so  that  the  reflection 
coefficients  can  be  found. 

In  the  next  step,  a  periodic  BBC  model  is  proposed,  from  which  we  can  derive  an 
expression  for  the  reflection  coefficients.  In  the  final  step,  the  optimal  parameters  for  the 
IBC  model  are  obtained  by  minimizing  the  mean  squared  error  between  the  two  sets  of 
reflection  coefficients  based  on  the  genetic  algorithm.  These  steps  are  described  in  detail 
below. 

A.  Periodic  IBC  model 

The  equivalent  BBC  relating  the  tangential  electric  and  magnetic  fields  for  a  planar 
coated  surface  can  be  written  as  [1]: 


We  shall  adopt  this  model  for  the  corrugated  problem  due  to  its  simplicity  and  usefulness 
for  our  applications.  The  model  parameters  will  then  be  optimized  to  emulate  the 


properties  of  the  exact  structure.  Note  that  since  the  corrugated  surface  exhibits 
anisotropic  scattering  characteristics,  the  equivalent  EBC  must  also  in  general  be 
anisotropic.  Therefore,  the  cross  impedance  terms  and  Zzz  are  kept  in  our  formulation 


to  assess  their  importance.  The  boundary  impedance  Zzz,  Zzx,  Zxz,  and  are  in  general 
functions  of  incident  angle  and  spatial  position.  For  the  IBC  model  to  be  useful  for 
subsequent  electromagnetic  analysis,  however,  it  is  much  preferable  to  model  the 
boundary  impedances  as  functions  of  spatial  position  only.  We  cannot  prove  theoretically 
the  existence  of  such  a  model  for  an  arbitrary  corrugated  structure.  Instead,  the 
applicability  and  limitation  of  this  approach  will  be  explored  numerically  in  Section  III. 

In  our  EBC  model,  the  periodic  grating  structure  with  period  p  as  shown  in  Fig.  2  (a) 
is  replaced  by  an  equivalent  impedance  boundary  condition  which  also  has  a  period  p,  as 
illustrated  in  Fig.  2  (b).  Each  surface  impedance  term  can  be  expanded  into  a  Fourier 
series.  Since  the  cross  impedance  terms  Zu  and  Z XT  are  usually  very  small,  we  shall  treat 
them  as  constants  and  only  expand  the  impedances  Z,x  and  Zxz  as 
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Therefore,  to  fully  describe  the  IBC  model,  we  must  determine  the  Fourier  series 
coefficients  {an}  and  {bn}. 


B.  Solution  to  the  forward  problem  of  scattering  by  the  IBC  model 

Next,  we  derive  the  reflection  coefficients  resulting  from  the  plane  wave  scattering 
from  the  IBC  model  given  above.  Under  plane  wave  incidence  where 

klx0  =  k0  sin 0 sirup  ,  k‘y0  =  -k0  cos#  and  k[0  =  k0  s'mOcostp  , 


each  component  of  the  tangential  electric  and  magnetic  fields  at  the  impedance  surface 
can  be  expanded  into  a  sum  of  Floquet  harmonics  [7].  For  example,  the  tangential  electric 
field  in  the  z  direction  is  expanded  as 

E  =  E‘e~jk,°x  +  V  Er  e~JkmX  (3) 
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where 

.  .  2k 

kxn  =kx0  +n — 

P 

is  the  propagation  constant  of  the  nth  order  harmonic  along  the  x-direction.  This  Floquet 
harmonic  is  a  reflected  wave  with  propagation  constants  (km,  k  yn,  k2),  where 
kr  =Jkn  —  k]  ~-~k^~  with  the  square  root  taken  as  positive  real  or  negative  imaginary. 
The  superscripts  i  and  r  denote  respectively  the  incident  and  reflected  field  throughout 
this  paper.  The  harmonic  term  ejc*e~jk‘z  is  suppressed  in  (3)  and  y  is  set  to  zero  at  the 
impedance  surface  for  convenience.  Assuming  that  the  coefficients  of  the  higher  order 
Floquet  harmonics  are  negligible,  and  applying  this  to  other  tangential  field  components, 
we  get 

N 

E  =  F'e~jk'oX  +  FTne~,kmX  (4) 

n=-N 

where  F  can  be  Ez,  Hz,  Ex  or  Hx,  and  N  is  a  positive  integer.  Substituting  (2)  and  (4)  into 
(1),  and  matching  the  coefficients  of  the  exponential  terms,  we  obtain  a  set  of  equations 
relating  E^,  Hxn  and  Hzn 

E'.Sn0  +  El  =  z=h;  +  +  a-H'‘  + 

m=-N 


(5) 


for  n  =  -N  to  N.  4io  is  the  Kronecker  delta.  This  set  of  equations  can  be  written  in  matrix 
form  as 
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Following  the  same  steps,  we  also  arrive  at  the  relationship  for  Exn,  Hxn  and  Hzn 
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For  a  plane  wave,  the  tangential  components  of  the  fields  in  the  x-direction  Ex,  Hx  can  be 
expressed  in  terms  of  Ez  and  Hz  as: 
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Combine  (6),  (7)  and  (8),  the  matrix  relationship  between  the  incident  and  reflected  fields 
is  written  as 
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In  (9),  the  z  components  of  the  scattered  field  can  be  directly  related  to  those  of  the 
incident  field.  We  can  therefore  define  the  reflection  coefficients  for  the  different 


polarizations  as: 


well  known  that  in  general  layered  coating  problems,  there  is  not  a  preferred  position  for 
the  impedance  surface.  The  solution  can  sometimes  be  improved  by  applying  the  IBC  at  a 
position  other  than  a  natural  interface  in  the  structure  [4],  Therefore,  by  including  the 
position  of  the  EBC  surface  as  an  additional  tuning  parameter  in  our  IBC  model,  we  can 
further  improve  the  accuracy  of  the  model. 


C.  Genetic  algorithm  to  determine  the  optimal  IBC  parameters 


In  the  genetic  algorithm,  the  parameters  to  be  optimized  are  first  encoded  into  binary 
form.  A  set  of  the  encoded  parameters  is  known  as  a  chromosome.  The  basic  idea  of  GA 
is  to  generate  a  pool  of  chromosomes,  discard  the  bad  ones,  keep  the  best  ones  and  let 
them  evolve  to  produce  better  chromosomes.  The  evaluation  of  each  chromosome  is 
performed  by  a  cost  function  which,  in  this  case,  is  chosen  to  be  the  mean  squared  error 
between  the  MoM  computed  reflection  coefficients  and  those  solved  using  (9)  and  (10) 
with  the  (ZK,  an,  bn,  Z„,  d)  parameters  decoded  from  the  corresponding  chromosome. 
Chromosomes  in  the  pool  are  ranked  according  to  the  cost  function.  The  best  ones  are 
selected  in  pairs  to  act  as  parents  of  the  next  generation.  Reproduction  of  children 
chromosomes  is  based  on  specific  rules  of  heredity  and  mutation.  The  process  of 
selection  and  reproduction  is  repeated  until  a  set  of  satisfactory  parameters  is  found  or  the 
generation  limit  is  reached.  The  flow  chart  of  the  genetic  algorithms  is  shown  in  Fig.  3. 
Detailed  discussion  of  the  genetic  algorithms  can  be  found  in  [8]. 

In  the  IBC  model  for  the  corrugated  coating,  the  parameters  to  be  optimized  are  the 
coefficients  of  the  Fourier  expansion  an  and  bn,  the  cross  impedances  Z„,  Zxx  and  the 
position  of  the  impedance  surface  d.  Each  of  the  parameters  a,„  bn,  Z~  and  ZTV  consists  of 
a  real  part  and  an  imaginary  part.  We  assume  a  symmetric  structure  so  that  an=a.n,  bn=b.n. 
For  an  approximation  truncated  to  the  Mh  order,  the  total  number  of  real  numbers  is 
4;V+9.  The  number  of  bits  contained  in  each  parameter  B  is  adjustable.  If  B  is  too  large, 
the  convergence  of  GA  will  be  slow.  If  B  is  too  small,  the  accuracy  of  the  calculation  will 
suffer.  In  the  examples  given  in  this  paper,  we  choose  B  =  8  to  be  efficient  in  both  speed 
and  accuracy.  In  order  to  encode  the  unknown  parameters  into  binary  form,  the  minimum 
and  maximum  possible  values  of  each  parameter  are  required.  For  example,  the  values 


for  the  real  and  imaginary  parts  of  the  Fourier  series  an  and  bn  are  estimated  to  be  in  the 
range  from  -3  r)0  to  3t]0,  which  is  found  to  be  reasonable  in  the  numerical  examples.  Thus 
the  8-bit  binary  00000000  denotes  -3  tj0  and  11111111  represents  3  r]0.  Z„  and  Z«  are 
relatively  small  and  their  real  and  imaginary  parts  vary  from  -0.1  rj0  to  0.1  T]0.  The 
distance  d  is  limited  between  the  upper  and  lower  boundaries  of  the  coating  so  that  the 
resulting  H3C  model  will  not  cause  any  ambiguity  in  its  applications. 

In  the  beginning  of  the  genetic  algorithm,  a  number  of  chromosomes  are  randomly 
generated.  Each  chromosome  is  decoded  into  parameters  Z7z,  an,  bn,  Zxx  and  d.  The 
reflection  coefficients  are  then  computed  using  (9)  and  (10).  The  cost  function  gives  the 
mean  squared  error  between  these  reflection  coefficients  and  their  corresponding  MoM 
solution: 

8.<t> 

where  R0  denotes  the  MoM  solution  of  the  reflection  coefficients  at  a  specific  observation 
angle  (0,  <p)  and  Pl5  P2  is  the  polarization  TEZ  /  TMZ  (or  V  /  H).  The  fitness  value  of  each 
chromosome  is  given  by 

f(Cl)  =  c[-c2Cost(C,) 

where  c\  and  c2  are  constants  and  C,  is  the  ith  chromosome  in  the  population.  This  fitness 
value  is  used  in  ranking  the  chromosomes  and  selecting  of  parents  for  the  next  generation 
[8-11].  There  are  several  standard  ways  of  selection.  In  this  paper,  the  roulette  wheel 
selection,  in  which  the  probability  of  each  chromosome  to  be  selected  is  proportional  to 


its  fitness  value,  is  used. 


After  two  chromosomes  are  selected,  they  mate  to  generate  children.  This  is  realized 
by  the  process  of  crossover,  in  which  a  break  point  is  randomly  chosen  in  the 
chromosomes  and  the  two  chromosomes  are  switched  at  that  point.  Mutation  is  imposed 
at  this  point  so  that  new  genes  appear  in  the  next  generation.  The  mutation  rate,  which  is 
the  portion  of  bits  to  be  randomly  changed,  is  also  an  important  parameter  in  GA. 
Experiments  show  that  a  mutation  rate  of  5-8%  is  often  efficient  in  the  calculation. 

The  process  of  evaluation,  selection  and  reproduction  is  repeated  until  a  desired  mean 
squared  error  is  achieved  or  a  maximum  generation  is  reached.  For  a  population  of  400 
chromosomes,  the  0th  order  EBC  (ie.  N  =  0)  takes  10-20  generations  to  converge  to  the 
optimum  while  the  2nd  order  EBC  takes  200-400  generations. 

III.  Numerical  Examples 

In  this  section,  some  examples  are  presented  to  demonstrate  the  effectiveness  of  the 
method.  The  first  example  is  a  deep,  triangular  grooved  grating  with  relatively  small 
period.  The  geometry  of  one  cell  of  the  grating  is  shown  in  Fig.  4  (a)  where  the  period  p 
=  0.067 Ao  and  hi  =  0.22  A#,  h2  =  0.017  A 0.  The  coating  material  is  MagRAM  with  material 
constants  er  =  14.35  -  y0.28  and  pr  =  1.525  -  ;'1.347  at  the  frequency  of  10  GHz.  17 
observation  angles  are  selected  which  include  normal  incidence  and  the  combination  of  9 
=  20°,  40°,  60°,  80°  and  0  =  0°,  30°,  60°,  90°.  In  this  example,  we  set  and  Z~  in  (1)  to 
zero  and  N  =  0  in  (4)  to  make  the  model  comparable  with  TEBC.  The  co-polarization 
reflection  coefficients  for  the  H-pol  and  V-pol  incidence  are  plotted  in  Figs.  5(a)  and 
5(b),  respectively,  and  the  H-V  cross-polarization  reflection  coefficients  are  plotted  in 
Fig.  5(c).  In  the  figures,  the  x  axis  is  divided  into  sections  of  different  incident  angle  9. 


which  varies  from  5  to  85  degrees  in  steps  of  10  degrees.  In  each  of  the  0  section,  the 
grating  angle  <p  varies  from  0  to  90  degrees.  The  matching  of  the  reflection  coefficients 
between  the  GA  approach  and  MoM  solution  is  good  at  most  incident  angles,  even  near 
grazing  incidence.  The  value  of  d  is  found  to  be  0.027/Jo  from  the  tip  of  the  groove,  or 
0.21^  above  the  ground  plane,.  The  TEBC  result  is  also  generated  by  using  the  reflection 
coefficient  at  normal  incidence  to  derive  the  equivalent  boundary  condition.  The 
impedance  surface  is  placed  at  the  plane  of  the  conductor  backing.  It  can  be  seen  that  the 
TIBC  results  deviate  significantly  from  the  reference  solution  away  from  normal 
incidence.  With  the  same  complexity  of  the  boundary  condition,  GA  achieves  a  much 
better  matching  because  more  observation  points  are  used  in  the  modeling  and  because  of 
the  additional  degree  of  freedom  in  the  position  of  the  EBC  surface. 

Next,  we  compare  the  IBC  approximations  with  and  without  the  cross  impedance 
terms  and  Z^.  The  structure  is  a  rectangular  groove  as  shown  in  Fig.  4  (b)  with  a 
period  p  =  QAlAo  and  a  groove  depth  of  hi  =  h2  =  0.042Ao.  The  material  constants  are  er  = 
8.3-y2.4  and  pr  =  2-;0.9  at  10  GHz.  The  same  observation  points  are  used  as  in  the 
previous  example.  The  0th  order  (N  =  0)  EBC  is  determined  and  the  comparison  between 
the  reflection  coefficients  is  illustrated  in  Fig.  6.  With  the  cross  impedance  included,  the 
accuracy  of  the  approximation  is  improved. 

In  the  third  example,  the  EBC  of  different  orders  are  obtained  for  the  structure  shown 
in  Fig.  4  (b)  where  p  =  0.42  Ao,  hi  =  0.25  Ao  and  h2  =  0.17  Ao.  The  coating  material  is  the 
same  as  that  in  example  2  but  the  period  is  much  larger  and  the  groove  is  deeper.  The 
reflection  coefficients  predicted  by  the  0th  order  and  2nd  order  model  are  plotted  in  Fig. 
7.  While  the  approximation  by  the  0th  order  model  is  fairly  satisfactory,  the  2nd  order 


model  further  improves  the  result  and  the  matching  is  better  at  most  incident  angles.  The 
price  of  the  improvement  is  the  computation  time.  For  the  Oth  order  modeling,  it  takes 
only  a  few  minutes  for  the  GA  to  converge  while  the  2nd  order  IBC  takes  more  than  1 
hour  on  an  SGI  02  workstation  (R10000/155MHz).  Another  consequence  as  the  order  of 
the  model  is  increased  is  that  the  resulting  IBC  will  show  more  spatial  variation.  This 
implies  that  when  the  IBC  model  is  utilized  in  subsequent  analysis  using  numerical 
electromagnetics  solvers,  the  impedance  surface  must  be  divided  more  finely  to  faithfully 
describe  the  IBC.  This  will  lead  to  a  higher  computation  cost.  Thus  the  higher  order 
model  is  not  recommended  unless  the  Oth  order  one  is  intolerable  or  the  period  is  large 
compared  to  the  wavelength.  Generally  speaking,  the  EBC  model  can  be  improved  by 
increasing  the  model  order,  whether  or  not  the  cross  impedance  terms  exist.  But  with  the 
cross  impedance  terms,  the  required  model  order  is  usually  smaller  than  that  without 
them. 

We  further  investigate  a  structure  with  a  period  larger  than  half  a  wavelength,  which 
results  in  higher  order  Floquet  mode  reflection  at  some  incident  angles.  The  rectangular 
groove  shown  in  Fig.  4(b)  has  a  period  p  =  0.852o  and  hi  =  /i?  =  0.42  Ao.  The  material 
constant  are  er  =  10.5-/2.2  and  jur  =  2-;0.3.  The  Oth  and  1st  order  reflection  coefficients 
are  plotted  in  Figs.  8(a)  and  8(b),  respectively.  It  is  shown  that  the  Floquet  modes  are  also 
well  characterized. 

In  the  final  example,  we  investigate  the  limitation  of  the  BBC  model.  We  consider  a 
triangular  groove  shown  in  Fig.  4(a)  with  p  =  0.33  Ac  and  /j/  =  /z2  =  0.0820.  The  optimal 
IBC  model  is  found  using  the  genetic  algorithm  for  different  coating  materials.  A  second 
order  IBC  model  with  cross  impedance  terms  is  used  and  the  optimal  model  parameters 


are  determined  by  running  GA  to  convergence.  After  the  model  is  found,  the  root  mean 
squared  (RMS)  error  of  the  EBC-predicted  reflection  coefficients  over  the  selected  angles 
are  computed.  Fig.  9  shows  the  RMS  error  of  the  optimal  IBC  model  as  a  function  of  er' 
and  er”,  which  are  the  real  and  imaginary  parts  of  the  coating  relative  permittivity  er.  /ur  is 
set  to  1  for  all  the  coatings.  We  observe  that  the  IBC  model  works  best  for  high  contrast, 
high  loss  materials.  For  low  contrast  or  low  loss  materials,  the  model  error  can  be  large. 
This  behavior  is  very  similar  to  conventional  IBC  models  for  planar  coatings. 

We  now  apply  our  derived  IBC  model  to  a  3D  scattering  problem.  Consider  the 
comer  reflector  as  shown  in  Fig.  10.  The  monostatic  radar  cross  section  (RCS)  is 
calculated  for  both  the  uncoated  reflector  and  that  coated  with  the  MagRAM  structure 
described  in  example  1.  The  groove  of  the  coating  is  either  parallel  or  perpendicular  to 
the  incident  direction.  Both  cases  are  computed  for  comparison.  Note  that  the  solution  for 
such  a  structure  is  very  complicated  if  we  try  to  use  the  exact  MoM  formulation.  Instead, 
we  use  the  0th  order  impedance  boundary  condition  obtained  from  example  1  to  replace 
the  corrugated  absorber.  We  assume  that  the  size  of  the  plate  remain  the  same  after  the 
EBC  replacement.  The  RCS  is  computed  using  FISC  [12],  which  is  a  3D  MoM  code 
based  on  the  fast  multipole  method  [13].  Comparisons  of  the  RCS  at  several  elevation 
angles  are  shown  in  Fig.  1 1  for  both  the  H-  and  V-polarizations.  The  result  shows  the 
effect  of  coating,  which  lowers  the  overall  RCS  level  for  both  polarizations.  We  further 
observe  that  a  30  dB  RCS  reduction  can  be  achieved  for  both  polarizations  over  the  range 
of  elevation  angles  from  20  to  75  degrees  if  the  grating  is  oriented  parallel  to  the  incident 


wave. 


IV.  Conclusions 


In  this  paper,  an  impedance  boundary  condition  model  is  derived  based  on  the  genetic 
algorithm  to  approximate  arbitrary  corrugated  coating  structures  in  scattering  problems. 
The  periodic  structure  is  replaced  by  a  periodic  IBC  on  a  virtual  surface.  The  boundary 
impedance  and  the  position  of  the  surface  are  optimized  by  matching  the  reflection 
coefficients  to  the  rigorous  numerical  solution  at  a  number  of  incident  angles.  Similar  to 
traditional  IBC  models,  this  approach  is  most  effective  when  the  coating  material  is  high- 
loss  and  of  high  contrast.  The  resulting  IBC  model  generated  by  this  algorithm  can  be 
incorporated  into  an  existing  computational  electromagnetics  code  to  assess  the 
performance  of  the  corrugated  coating  when  mounted  on  complex  platforms. 

Compared  with  other  IBC  approaches,  the  method  described  above  has  several 
advantages.  First,  the  boundary  impedance  is  assumed  to  be  anisotropic  so  that  the  same 
model  can  be  applied  to  oblique  incidence  from  any  arbitrary  angles.  Second,  it  is 
possible  to  build  in  spatial  variation  of  the  boundary  impedance  by  adjusting  the  number 
of  terms  used  in  the  Fourier  series  expansion.  By  using  more  terms,  the  IBC  model  can  be 
made  more  accurate.  In  addition,  the  position  of  the  impedance  surface  can  also  be 
optimized.  By  solving  for  the  best  position  of  the  impedance  surface  as  one  of  the  model 
parameters,  the  accuracy  of  the  model  can  be  improved. 

Numerical  experiments  show  that  the  IBC  approximation  can  be  improved  if  some  of 
the  parameters  of  the  genetic  algorithm  are  carefully  chosen.  These  parameters  include 
the  incident  angles  at  which  the  rigorous  solution  is  obtained,  the  range  of  each  model 
parameter  and  the  mutation  rate,  etc.  The  genetic  algorithm  can  also  be  accelerated  with 
carefully  chosen  parameters  and  a  well  designed  cost  function. 
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Fig.  1.  Scattering  from  a  corrugated  coating  structure  backed  by  a  conducting  surface 


(a)  Original  structure 


(b)  EBC  Approximation 


Fig.  2.  Equivalent  impedance  boundary  problem 


Incident  angle 
(c)  H-V  cross-polarization 

Fig.  5.  Comparison  of  reflection  coefficients  versus  angle  among  (i)  the  exact 
MOM  result,  (ii)  the  EBC  derived  from  GA  and  (iii)  the  TIBC  result 
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Fig.  6.  Effect  of  incorporating  the  cross  impedance  terms  in  the  IBC  model  generated  from  GA 
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Fig.  7  (d)  Phase,  V-polarization 
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Fig.  8  Higher  order  reflection  coefficients  simulated  by  the  EBC  model 


Fig.  10.  Geometry  of  the  comer  reflector 
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ABSTRACT 

A  radar  cross  section  (RCS)  interpolation  technique  in  both  frequency  and 
aspect  is  proposed  for  the  efficient  prediction  of  radar  signatures  from 
computational  electromagnetics  data.  Our  approach  is  based  on  a  multiple-arrival 
model  for  the  induced  current  on  the  target.  The  model  parameters  are 
determined  by  an  adaptive  feature  extraction  (AFE)  algorithm,  which  uses  an 
iterative  search-and-extract  procedure  to  find  the  individual  model  features. 
Random  frequency  and  aspect  sampling  is  used  to  circumvent  the  ambiguity  in 
selecting  the  features.  Numerical  examples  are  presented  to  test  the  interpolation 
algorithm.  It  is  found  that  sufficient  accuracy  in  the  predicted  radar  features  can 
be  achieved  even  when  the  original  computed  data  is  sampled  at  5:1  below  the 
Nyquist  criterion  in  either  frequency  or  aspect.  The  algorithm  is  also  applied  to 
efficiently  predict  the  radar  image  of  the  benchmark  VFY218  airplane  at  UHF 
band  with  good  results. 


I. 


Introduction 


Accurate  prediction  of  the  radar  scattering  from  complex  targets  requires 
solving  the  electromagnetic  boundary-value  problem  using  numerical  methods. 
This  is  a  computationally  intensive  task  when  targets  of  large  electrical  sizes  are 
involved.  Recently,  fast  iterative  solvers  have  emerged  based  on  the  fast  multipole 
method  [1,2].  These  solvers  have  much  lower  computational  complexity  than 
traditional  approaches.  However,  for  frequency-aspect  radar  cross  section  (RCS) 
calculations,  the  solver  has  to  be  executed  repeatedly  for  each  angle  and  each 
frequency.  To  generate  RCS  data  with  as  few  electromagnetic  computation  points 
as  possible,  a  number  of  model-based  interpolation  or  extrapolation  approaches 
have  been  developed  to  date  [3-10].  In  this  work,  we  address  the  problem  of  RCS 
interpolation  in  frequency  and  aspect  based  on  a  multiple-arrival  model  in 
conjunction  with  an  adaptive  feature  extraction  (AFE)  algorithm. 

Model-based  techniques  are  based  on  the  assumption  that  all  the  unknown 
samples  to  be  interpolated  obey  the  same  reduced-order  model  as  the  known 
samples,  from  which  the  model  coefficients  can  be  determined.  The  key  in  the 
success  of  a  model-based  interpolation  algorithm  is  a  good  model  of  the  physical 
observable  to  be  interpolated.  The  most  commonly  used  model  is  the  rational 
function  model,  which  is  applicable  in  the  low  frequency  region.  Here,  we 
consider  a  multiple-arrival  model  for  the  induced  current.  This  model  is 
motivated  by  the  well-known  behaviors  of  scattering  mechanisms  at  high 
frequencies  [11],  and  has  been  applied  previously  to  the  extrapolation  of  RCS  in 
the  frequency  and  angular  domains  with  good  success  [8-10]. 


To  obtain  the  model  parameters  from  the  computed  data,  we  adopt  an 
algorithm  termed  adaptive  feature  extraction  (AFE).  AFE  can  be  considered  as  a 
generalization  of  the  CLEAN  algorithm  [12],  and  is  similar  to  the  adaptive  joint 
time-frequency  technique  [13-15]  and  the  matching  pursuit  algorithm  [16]  in  the 
iterative  manner  it  performs  the  parameterization.  We  have  previously  applied  it 
to  construct  the  inverse  synthetic  aperture  radar  (ISAR)  image  from  radar 
measurement  data  that  was  undersampled  in  the  aspect  dimension  [17].  The 
essential  idea  of  the  AFE  algorithm  is  to  search  and  extract  out  individual 
scattering  features  from  the  data  set  one  at  a  time.  During  each  iteration  the 
strongest  feature  is  identified  and  removed  from  the  original  data.  The  procedure 
is  then  iterated  until  the  data  is  well  parameterized  by  the  feature  set.  In  this 
manner,  the  interference  between  different  scattering  features,  which  is  significant 
for  undersampled  data,  can  be  largely  avoided.  Since  the  features  in  our  model 
are  exponential  functions  of  frequency  and  aspect,  we  use  random  sampling 
during  the  collection  of  the  original  computation  data  in  order  to  avoid  the 
ambiguity  in  selecting  the  strongest  feature.  This  approach  is  similar  to  the 
random  array  concept  that  uses  highly  thinned  but  randomly  spaced  elements  to 
avoid  grating  lobes  [18]. 

This  paper  is  organized  as  follows.  In  Sec.  2,  we  introduce  the  multiple- 
arrival  model  for  the  current  and  the  AFE  algorithm  to  develop  a  one-dimensional 
(1-D)  frequency  interpolation  scheme.  A  numerical  example  is  given  to  verify  our 
approach  and  show  that  the  Nyquist  sampling  criterion  can  be  overcome  with  the 
extra  information  provided  by  the  model-based  approach.  In  Sec.  3,  the 


interpolation  algorithm  is  extended  for  2-D  frequency-aspect  interpolation.  Some 
numerical  examples  are  given  to  illustrate  the  performance  of  the  algorithm.  The 
algorithm  is  also  applied  to  efficiently  predict  the  ISAR  image  of  the  benchmark 
VFY218  airplane  at  UHF  band.  Finally,  some  discussions  and  conclusions  on  the 
interpolation  algorithm  are  given  Sec.  4. 

II.  1-D  Frequency  Interpolation 

2.7.  Multiple- Arrival  Model  and  AFE  Interpolation  Algorithm.  Before  presenting 
our  interpolation  algorithm,  we  shall  briefly  discuss  the  sampling  criterion  for 
frequency  response  in  electromagnetic  scattering  problems.  Nyquist  sampling 
theorem  states  that  if  the  sampling  rate  is  higher  than  the  highest  Fourier 
component  of  a  signal,  the  signal  can  be  perfectly  reconstructed  by  interpolation 
with  the  sine  function.  In  the  context  of  electromagnetic  scattering,  this  means 
that  the  frequency  sampling  interval  for  the  scattered  far  field  should  be  less  than 
1  /(W  )>  where  Tmax  corresponds  to  the  total  extent  of  the  time-domain  response. 
Strictly  speaking,  due  to  the  finite  frequency  bandwidth  of  the  data,  the  actual 
time  extent  of  the  signal  does  not  have  finite  support  and  a  perfect  reconstruction 
is  not  guaranteed.  However,  an  approximate  Nyquist  criterion  can  usually  be 
defined  in  practice  for  the  sampling  density  in  frequency. 

To  overcome  the  Nyquist  sampling  limitation,  extra  information  must  be 
provided  by  an  a  priori  model  based  on  the  scattering  physics.  We  utilize  a 
multiple  time-of-arrival  model  of  the  induced  current: 

W=^BP  expf  -  jlrftp) 


(1) 


The  essential  idea  of  this  model  is  illustrated  in  Fig.l.  It  assumes  that  the  current 
on  the  target  surface  is  induced  by  various  incident  wave  mechanisms  that  arise 
from  the  direct  incident  field  as  well  as  the  multiply  scattered  fields  from  other 
parts  of  the  target.  Each  incident  mechanism  has  amplitude  Bp,  which  is  assumed 
to  be  frequency  independent,  and  a  time-of-arrival  tp.  Thus  each  basis  function  in 
(1)  is  an  exponential  function  with  linear  phase  in  frequency.  This  is  what  we 
shall  also  refer  to  as  a  “feature”  in  the  data.  This  model  has  been  found  to  work 
well  in  modeling  the  scattering  mechanisms  of  complex  targets  [8-10].  It  should 
be  noted  that  the  induced  current  instead  of  the  scattered  far  field  is  chosen  for  the 
interpolation.  This  is  because  the  scattered  field  in  general  requires  a  much  higher 
model  order  to  parameterize.  Furthermore,  each  current  element  usually  contains 
more  localized  interactions  in  time  than  the  total  scattered  field  and  therefore  has 
a  looser  Nyquist  sampling  criterion.  The  price  we  pay  is  that  the  interpolation 
needs  to  be  performed  for  all  the  current  elements  on  the  target. 

If  the  original  frequency  sampling  is  dense  enough,  the  basis  functions  in  (1) 
are  orthogonal  and  the  model  parameters  Bp  and  tp  can  be  determined  using  the 
Fourier  transform.  However,  when  the  data  is  undersampled,  the  basis  functions 
are  no  longer  orthogonal.  To  avoid  the  interference  between  the  non-orthogonal 
bases,  we  use  the  iterative  AFE  procedure  instead  of  the  Fourier  transform  to 
carry  out  the  parameterization.  To  find  parameters  Bp  and  tp,  we  first  project  the 
sampled  frequency  response  onto  the  complex  conjugate  of  the  model  bases  for 
all  possible  values  of  t.  We  then  select  as  the  best  basis  the  one  with  the 


maximum  projection  value  Bp  and  the  corresponding  time-of-arrival  parameter  tp. 
This  is  denoted  by  the  following  equation: 


Bp  =max(yp(/),exp0'2^)) 


(2) 


where  the  subscript  p  denotes  it  is  at  the  pth  stage  of  the  iterative  procedure  and 
the  inner  product  is  defined  as: 


Once  the  strongest  feature  is  found,  a  remainder  signal  is  generated  by  subtracting 
out  that  feature: 


(/,  )  =  Jp(fi)~Bp  exp(-j2ftf.tp )  (4) 

The  above  process  is  iterated  to  extract  out  as  many  features  as  needed  to 
sufficiently  parameterize  the  signal,  i.e.,  until  the  remainder  signal  reaches  a 
preset  threshold.  Note  that  since  the  stronger  features  are  extracted  out  before  we 
begin  the  search  for  the  weaker  features  in  the  iteration  process,  the  interference 
from  the  stronger  features  on  the  weaker  features  is  reduced  significantly. 

An  additional  problem  exists  due  to  the  undersampling  of  the  exponential 
bases  used  in  our  model.  As  an  example,  if  J(f)  consists  of  a  single  exponential 
function  exp(-j2nfto),  we  expect  the  inner  product  in  (2)  to  give  the  strongest  peak 
at  tp  =  t0  with  amplitude  Bp  =  1 .  However,  if  the  frequency  sampling  is  uniform 
with  sampling  interval  $f,  the  strongest  peak  will  occur  not  only  at  tp  =  to  but  also 
at  tp  =  t0  +2n7t/Sf.  These  periodic  lobes  cause  ambiguity  in  determining  the  true 
time-of-arrivals.  To  overcome  this  problem,  we  utilize  random  sampling  to 
weaken  the  repeated  lobes.  The  idea  is  similar  to  random  arrays  where  random 


element  spacing  is  used  to  circumvent  the  antenna  grating  lobe  problem  in 
designing  highly  thinned  arrays  [18].  The  random  sampling  in  frequency  can  be 
easily  realized  in  the  present  problem  without  computational  penalty.  We  choose 
the  frequency  sampling  points  fj  based  on  a  uniform  probability  distribution  in  the 
frequency  band  of  interest.  Then  once  the  current  is  parameterized  at  each  point 
on  the  scatterer  using  AFE,  the  interpolated  current  function  at  any  frequency 
within  the  band  can  be  calculated  using  equation  (1).  Finally,  the  far  field  at  the 
denser  frequency  sampling  can  be  obtained  by  integrating  the  induced  current 
over  the  target. 

2.2.  Numerical  Example.  To  verify  the  above  frequency  interpolation  scheme, 
we  consider  the  scattering  from  a  two-dimensional  structure  consisting  of  three 
circular  cylinders  as  an  example.  The  structure  is  shown  in  Fig.  2.  The  reference 
RCS  is  computed  at  71  points  from  0.3  GHz  to  0.65  GHz  using  the  method  of 
moments.  The  result  is  shown  as  the  solid  line  in  Fig.  3(a).  The  spiky  behavior  of 
the  frequency  response  indicates  that  the  scattered  field  is  sampled  very  close  to 
the  Nyquist  rate.  For  comparison,  we  use  the  values  at  18  equally  spaced  points 
and  carry  out  a  current-based  interpolation  using  simple  spline  fitting.  The  result 
is  shown  as  the  dash-dotted  line  in  Fig.  3(a).  Deviations  between  the  two  results 
can  be  seen.  Next,  we  use  the  AFE  algorithm  to  carry  out  the  interpolation. 
Instead  of  18  equally  spaced  points,  we  randomly  select  18  points  from  the 
original  71  points.  The  interpolated  result  is  plotted  as  the  dashed  line  in  Fig. 
3(b),  which  agrees  with  the  reference  calculation  much  better  than  the  current- 
domain  spline  interpolation.  This  contrast  is  further  demonstrated  when  we 


Fourier  transform  the  complex  frequency  responses  into  the  range  profiles 
displayed  in  Fig.  4.  In  the  reference  range  profile,  the  first  two  strong  peaks  in 
range  are  the  direct  scattering  from  the  two  cylinders  on  the  left.  Those  smaller 
peaks  further  down  range  represent  the  multiple  scattering  due  to  the  interactions 
among  the  cylinders.  We  can  see  that  most  of  the  features  from  the  AFE 
interpolation  coincide  with  the  reference  result  in  Fig.  4(b),  while  the  current- 
domain  spline  interpolation  gives  strong  artifacts  among  the  real  scattering 
features  in  Fig.  4(a).  The  only  peaks  correctly  predicted  by  the  current  spline  are 
the  two  direct  reflection  peaks  off  the  two  left  cylinders.  This  is  because  the 
current  elements  contributing  to  those  two  peaks  contain  only  the  physical  optics 
component.  Since  each  current  element  has  only  a  single  time-of-arrival,  the 
Nyquist  sampling  criterion  is  very  loose  for  these  elements.  This  is  why  the 
current  at  these  positions  can  be  easily  interpolated  even  with  the  spline  function. 
Over  the  rest  of  the  target,  this  is  clearly  not  the  case. 

III.  2-D  Frequency- Aspect  Interpolation 

3.  L  2-D  AFE  Interpolation  Algorithm,  We  now  extend  the  above  1-D  frequency 
interpolation  algorithm  to  2-D  frequency-aspect  interpolation.  Note  that  for 
iterative  electromagnetic  solvers,  the  solver  has  to  be  executed  for  each  new 
aspect  angle,  just  as  in  the  frequency  dimension.  We  again  utilize  the  multiple- 
arrival  model  for  the  induced  current  shown  in  Fig.  1 ,  except  we  consider  how  the 
current  varies  as  a  function  of  both  the  frequency  and  the  direction  of  the  incident 
wave.  We  shall  denote  the  down  range  direction  with  respect  to  the  incident  wave 


as  x  and  the  direction  perpendicular  to  x,  or  the  cross  range  direction,  as  y.  The 
first  scattering  point  on  the  target  due  to  the  pth  scattering  mechanism  is 
represented  by  (xp,  yp).  lp  represents  the  path  length  from  the  first  scattering  point 
to  the  observation  point  S.  It  is  zero  for  the  direct  incident  wave.  Next,  we  make 
a  key  assumption  about  the  induced  current  at  S  as  the  angle  of  the  incident  wave 
is  varied  from  the  nominal  direction  (solid  arrows)  by  an  amount  6  shown  by  the 
dashed  arrows  in  Fig.  1.  We  assume  that  the  intermediate  interaction  path  /p  for 
the  pth  scattering  mechanism  remains  unchanged  as  the  incident  angle  is  varied. 
This  approximation  was  first  introduced  for  ray  optical  fields  [19],  and  has  been 
applied  in  model-based  angular  extrapolation  by  us  previously  with  good  success 
[9].  Using  this  assumption,  we  arrive  at  a  2-D  version  of  the  multiple-arrival 
model  for  the  current  at  a  point  S  on  the  target: 

J(f .  0=2  BP  exp  cos  6  +  y p  sin  6  +  lp )] 

P  C 

(5) 

Under  further  assumption  of  small  observation  angles,  (5)  can  be  approximately 
written  as: 

J(f,Q)~'LBP  exp  [-  j^-{  ( xp  +lp)cosd  +  ypsin  6)] 

P  C  (6) 

=  'LBP  exp  [  -  j—(rp  cos  0  +  yp  sin  6)] 

P  c 

Note  that  (6)  reduces  to  (1)  for  (9=  0,  rp=  ctp,  as  expected.  In  contrast  to  (1),  the 
basis  functions  in  (6)  are  two-dimensional  exponential  functions  of  frequency  and 
angle.  Therefore  we  modify  the  AFE  algorithm  in  Sec.  2  by  replacing  (2)  by: 


-  (7) 


Bp  =  max  (jp(f,0),exp[j^-(rcos0  +  ysm0)]\ 

r=rp,y=yp\  y  C  / 

By  searching  out  the  parameters  Bp,  rp  and  yp,  the  AFE  technique  is  thus  extended 
for  2-D  interpolation  in  both  the  frequency  and  angular  domains. 

To  summarize,  our  goal  for  the  2-D  interpolation  algorithm  is  to 
efficiently  generate  dense  frequency-aspect  RCS  data  from  a  small  number  of 
data  points  calculated  using  a  fast  iterative  solver.  The  steps  are  first  to  select  a 
number  of  randomly  distributed  sample  points  within  the  frequency  and  angular 
band  of  interest  and  solve  the  scattering  problem  for  those  points.  Second,  the 
induced  current  is  interpolated  to  a  denser  sampling  grid  based  on  the  reduced- 
order  model  generated  from  APE  processing.  Third,  the  current  is  integrated  to 
generate  the  scattered  far  field  and  the  RCS.  In  the  following  numerical  tests,  we 
go  one  step  further  by  generating  the  2-D  ISAR  image  of  the  target  from  the 
interpolated  frequency-aspect  data  using  standard  ISAR  processing.  This  is  done 
to  facilitate  better  physical  interpretation  of  the  results. 

32.  Numerical  Example.  To  test  the  frequency-aspect  interpolation  algorithm, 
we  consider  a  2-D  cylinder-plate  structure  shown  in  Fig.  5.  The  diameter  of  the 
cylinder  is  4.2  m  and  the  length  of  the  plate  is  20  m.  The  origin  of  the  cylinder 
and  the  center  of  the  plate  are  separated  by  6.2  m.  The  frequency  band  of  interest 
is  from  0.3  GHz  to  0.65  GHz.  Fig.  6  shows  the  reference  ISAR  image  generated 
from  7 lx  81=5751  computed  points  in  the  frequency- aspect  plane.  The  image  has 
a  dynamic  range  of  40  dB.  In  this  structure,  strong  multiple  scattering 
mechanisms  between  the  cylinder  and  the  plate  dominate  the  backscattering.  We 
can  clearly  see  in  the  image  the  features  corresponding  to  the  direct  scattering 


from  the  cylinder  and  the  edge  points  of  the  plate,  as  well  as  the  other  range- 
delayed  features  corresponding  to  the  multiple  scattering  mechanisms.  (The 
detailed  description  of  the  different  mechanisms  can  be  found  in  [19].)  For 
comparison,  we  collect  15x15=225  equally  spaced  sampling  points  in  frequency- 
aspect  and  use  2-D  splines  to  interpolate  the  current.  Fig.  7  plots  the  resulting 
ISAR  image  from  the  spline-interpolated  frequency-aspect  data.  In  this  image, 
only  the  direct  scattering  features  (the  cylinder  and  the  edge  points  of  the  plate) 
are  correctly  predicted  in  position  while  all  the  higher-order  scattering  features  are 
completely  wrong. 

Next  we  construct  the  ISAR  image  using  the  AFE  algorithm.  We  choose  50 
points  randomly  from  the  original  5751  points  and  use  AFE  to  interpolate  the 
frequency-aspect  data  to  the  dense  sampling  grid.  The  resulting  ISAR  image  is 
plotted  in  Fig.  8(a).  We  repeat  this  process  using  AFE  interpolation  for  100,  200 
and  400  input  points,  and  the  resulting  images  are  shown  in  Figs.  8(b),  8(c)  and 
8(d),  respectively.  From  this  series  of  images,  we  can  see  that  the  strong 
scattering  features  are  quite  stable  throughout  while  the  weak  features  begin  to 
converge  as  the  number  of  input  points  is  increased.  Nearly  all  of  the  features  of 
the  interpolated  image  agree  well  with  those  in  the  reference  image  from  the 
brute-force  calculation.  We  do  notice  that  the  very  weak  features  are 
underpredicted  in  the  interpolated  result.  This  is  mainly  caused  by  the 
termination  threshold  in  the  parameterization  process  and  some  signal  energy  is 
lost.  In  this  example,  the  parameterization  is  stopped  when  the  maximum 
magnitude  of  the  remainder  signal  is  less  than  10%  of  the  incident  magnetic  field 


strength.  Another  possible  cause  for  the  noise  in  the  low  dynamical  ranges  is  that 
the  AFE  algorithm  itself  leads  to  biased  estimates  on  the  position  and  amplitude 
of  the  scattering  features.  This  bias  becomes  more  severe  for  the  weaker  features. 

Plotted  in  Fig.  9  is  the  correlation  index  between  the  images  generated  using 
AFE  and  the  reference  image  versus  different  number  of  input  points  used  in  the 
interpolation.  To  properly  reflect  the  contribution  of  the  weaker  features  to  the 
correlation  index,  we  use  the  dB-scaled  images  in  computing  the  correlation  index 
(a  constant  offset  is  added  to  the  images  to  avoid  negative  numbers  in  decibel). 
Fig.  9  shows  that  the  200-point  mark  is  the  turning  point  of  the  curve,  after  which 
the  image  quality  does  not  improve  further.  This  is  consistent  with  the  visual 
sensation  in  Fig.  8.  From  this  example,  we  believe  a  good  rule-of-thumb  limit 
that  can  be  handled  by  the  AFE  interpolation  algorithm  is  about  5:1 
undersampling  in  both  the  frequency  and  aspect  dimensions.  This  criterion  is 
established  based  on  a  very  challenging  target  containing  strong  multiple 
scattering. 

3.3.  Application  of_AFE_  Interpolation  to  VFY218  Benchmark.  The  interpolation 
algorithm  is  next  applied  to  predict  the  ISAR  image  of  the  benchmark  VFY218 
airplane  [20].  The  fuselage  length  of  the  airplane  is  15.33  meters  and  the 
maximum  width  measured  from  the  wing  tips  is  8.90  meters.  To  generate  its 
ISAR  image  at  a  center  frequency  of  400  MHz  with  bandwidth  267  MHz,  we  use 
the  fast  solver  FISC  [21]  on  a  Pentium  II  450MHz  computer.  The  total  number  of 
sampling  points  must  be  at  least  40  frequency  points  by  40  aspect  points  within 
the  40-degree  aperture  to  satisfy  the  Nyquist  sampling  criterion.  The  resulting 


range  and  cross  range  resolution  is  about  half  a  meter.  Since  the  calculation  for 
one  single  frequency-aspect  point  takes  about  3  hours  (with  about  80,000 
unknowns  in  the  moment  equation),  the  total  computation  time  would  be  200  days 
if  we  use  brute-force  calculation  to  generate  the  data.  Based  on  the  5:1 
undersampling  criterion  derived  in  the  previous  example,  we  only  compute  the 
scattering  problem  at  62  randomly  sampled  points  and  use  the  AFE  interpolation 
scheme  to  interpolate  the  data  to  all  40x40=1600  points.  The  computation  time  is 
about  8  days.  The  resulting  ISAR  image  at  the  130-degree  (from  nose-on)  look 
angle  is  plotted  in  Fig.  10.  For  comparison,  Fig.  11  shows  the  ISAR  image 
constructed  from  the  chamber  measurement  data  for  the  same  look  angle.  The 
target  outline  is  overlaid  on  the  measurement  image.  The  measurement  was 
carried  out  at  the  US  Navy  China  Lake  facility  on  a  1:30  scaled  model  at  the 
frequency  band  of  8  to  16  GHz  [20].  Comparing  Fig.  10  to  Fig.  1 1,  we  find  that 
all  the  features  in  the  measurement  image  are  well  predicted  in  the  simulated 
image  from  using  FISC  and  interpolation.  The  correlation  index  (on  a  dB  scale) 
is  88.2  %  between  the  two  images.  The  sources  of  discrepancy  in  this  case 
include  measurement  error,  computation  error  from  FISC,  differences  between  the 
actual  model  used  in  the  measurement  and  the  CAD  model  used  in  the  FISC 
calculation,  and  error  from  the  interpolation  algorithm.  There  is  also  an  additional 
error  in  the  measured  image  due  to  the  ISAR  formation  procedure  [22].  For  the 
wide  40-degree  observation  aperture  in  the  measurement,  the  scattered  field  data 
is  collected  on  a  polar  grid  in  k-space  and  polar  reformatting  must  be  performed 
to  generate  a  focused  ISAR  image  with  FFT.  This  process  introduces  additional 


error  and  results  in  a  more  diffused  resolution  cell.  The  predicted  image,  on  the 
other  hand,  is  generated  directly  on  a  rectangular  grid  with  rectangular  boundary 
in  k-space  using  the  parameterized  current  model  and  completely  circumvents  the 
polar  reformatting  step. 

As  a  second  case  study,  the  predicted  image  (generated  using  62  FISC- 
computed  points  and  AFE  interpolation)  for  another  look  angle  at  30  degrees  from 
nose-on  is  shown  in  Fig.  12.  Two  measurement  images  are  shown  in  Figs.  13  and 
14.  Fig.  13  is  from  the  measured  data  when  the  airplane  engine  inlet  ducts  are 
untreated  and  Fig.  14  is  from  the  measured  data  when  absorber  material  is 
inserted  into  the  inlet.  Since  the  CAD  model  used  to  carry  out  the  simulation  has 
sealed  engine  inlets,  the  cloud  over  the  right  wing  in  Fig.  1 3  due  to  the  inlet  return 
is  not  predicted  in  Fig.  12.  However,  the  other  features  agree  quite  well.  The 
more  diffused  spot  sizes  in  the  measured  image  are  due  to  the  polar-reformatting 
operation,  as  discussed  earlier.  Fig.  12  also  agrees  well  with  Fig.  14,  although 
the  inlet  contribution  is  not  completely  removed  from  the  measurement  in  the 
latter  case.  Furthermore,  there  appears  to  be  more  measurement  noise  in  Fig.  14. 

IV.  Summary  and  Discussions 

A  current-domain  AFE  interpolation  technique  has  been  developed  to 
generate  densely  sampled  frequency-aspect  RCS  data  from  a  sparse  set  of 
computed  data.  The  induced  current  on  the  target  is  modeled  by  a  multiple-arrival 
model  and  the  AFE  procedure  is  applied  to  extract  the  model  parameters. 
Random  sampling  of  the  original  data  is  used  to  avoid  the  ambiguity  problem 


associated  with  the  exponential  model  basis.  Numerical  results  have  been 
generated  to  test  this  approach.  Lastly,  the  algorithm  has  been  applied  in 
conjunction  with  a  fast  electromagnetic  solver  to  predict  the  high-fidelity  ISAR 
image  of  the  benchmark  VFY218  airplane  at  UHF  band.  The  resulting  images 
agree  well  with  chamber  measurement  data. 

The  AFE  approach  is  found  to  be  stable  and  robust.  Sufficient  accuracy  in  the 
predicted  image  features  can  be  achieved  even  when  the  original  computed  data  is 
sampled  at  5:1  below  the  Nyquist  criterion  in  either  frequency  or  aspect.  This 
limit  was  established  based  on  a  target  containing  very  strong  multiple  scattering. 
Therefore,  the  expected  time  savings  in  using  this  approach  is  about  25:1  in 
comparison  to  the  brute-force  computation.  The  AFE  algorithm  does  involve 
exhaustive  search  and  must  be  carried  out  for  every  current  element  on  the  target. 
However,  the  time  consumed  in  the  interpolation  is  still  relatively  insignificant 
when  compared  to  the  electromagnetic  computation  time. 

The  possible  sources  of  error  in  the  overall  interpolation  procedure  include 
imperfection  of  the  multiple-arrival  model  for  the  current  and  errors  in  the  AFE 
procedure.  The  former  includes  our  assumption  about  the  frequency  independence 
of  the  model  coefficients  as  well  as  the  small-angle  approximations  used  to  derive 
the  model  in  the  angular  dimension.  Errors  in  the  AFE  procedure  can  occur  when 
there  are  too  many  features  in  the  data  or  when  the  undersampling  is  too  severe. 
The  interference  from  the  weaker  features  affects  the  determination  of  the  strong 
features  in  the  projection  process.  These  biased  estimates  of  the  strong  features, 
when  propagated  to  the  later  stage  of  the  iteration  process,  prevent  the  accurate 


extraction  of  the  very  weak  features.  This  error  source  limits  the  dynamic  range 
of  the  final  interpolated  data. 
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Fig.2  Geometry  of  the  three-cylinder  target. 
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Fig.9  Image  correlation  index  between  the  interpolated  result  and  the 

reference  result  versus  the  number  of  points  used  for  interpolation. 


Fig.  10 ISAR  image  of  VFY-218  at  130  degrees  from  nose-on 
generated  from  interpolated  result  using  AFE  with  62 
FISC-computed  points. 
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Fig.  11  ISAR  image  of  VFY-218  at  130  degrees  from  nose-on 
generated  from  chamber  measurement  data. 
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Fig.  12 ISAR  image  of  VFY-218  at  30  degrees  from  nose-on 
generated  from  interpolated  results  using  AFE  with  62 
FISC-computed  points.  The  CAD  model  used  in  the 
computation  has  sealed  engine  inlets. 
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Fig.  1 3  ISAR  image  of  the  VFY218  at  30  degrees  Fig.14  ISAR  image  of  VFY-218  at  30  degrees 
from  nose-on  generated  from  chamber  from  nose-on  generated  from  chamber 

measurement  data.  The  measured  model  measurement  data.  The  measured  model 

has  open,  untreated  engine  inlets.  has  absorber  material  placed  in  the  inlets. 
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Abstract 

A  frequency-aspect  extrapolation  algorithm  is  proposed  to  accelerate  ISAR  image 
simulation  using  fast  multipole  solvers.  A  2-D  multiple-arrival  model  based  on  high- 
frequency  physics  is  proposed  to  parameterize  the  induced  currents  on  the  target.  A  2-D 
ESPRIT  algorithm  is  developed  to  estimate  the  model  parameters  from  a  limited  number 
of  computed  data  samples  in  frequency  and  aspect.  The  model  is  then  extrapolated  to 
other  frequencies  and  aspects  to  arrive  at  broadband,  wide-angle  RCS  data  for  ISAR 
image  construction.  This  algorithm  is  tested  using  a  canonical  cylinder-plate  structure  to 
evaluate  its  performance.  The  ISAR  image  of  the  benchmark  VFY-218  airplane  at  UHF 
band  is  then  predicted  using  the  fast  multipole  solver  FISC  and  the  2-D  extrapolation 
algorithm.  The  resulting  image  compares  favorably  with  that  obtained  from  chamber 
measurement  data. 

Key  words:  ISAR  Image  Simulation,  Frequency-Aspect  Extrapolation.  2-D  ESPRIT, 
Fast  Multipole  Solver 


I.  INTRODUCTION 


Inverse  synthetic  aperture  radar  (ISAR)  imaging  is  an  important  tool  for  radar 
signature  diagnostic  and  target  identification  [1,2].  In  order  to  simulate  the  ISAR  image 
of  a  target,  it  is  necessary  to  solve  the  electromagnetic  scattering  problem  at  multiple 
frequencies  and  angles.  Then  by  performing  a  two-dimensional  (2-D)  Fourier  transform 
on  the  resulting  frequency-aspect  radar  cross  section  (RCS)  data,  a  2-D  ISAR  image  of 
the  target  can  be  constructed.  In  this  work,  we  consider  ISAR  image  simulation  of 
complex  targets  using  a  full-wave  numerical  technique.  The  focus  of  our  attention  is  the 
class  of  iterative  solvers  based  on  the  fast  multipole  method  [3,4].  These  solvers  have 
much  lower  computational  complexity  than  traditional  moment  method  for  scattering 
problems  at  a  single  frequency  and  a  single  observation  angle.  For  multiple  frequency- 
aspect  RCS  calculations,  however,  the  solver  has  to  be  executed  repeatedly  for  each  angle 
and  frequency. 

To  reduce  the  computation  burden  incurred  by  multiple  frequency-aspect 
calculations,  the  concept  of  model-based  parameter  estimation  can  be  applied  to  populate 
the  required  data  set  from  a  sparse  set  of  computed  data  [5-12].  In  particular,  we  have 
recently  developed  a  1-D  frequency  extrapolation  algorithm  [10.12]  and  a  1-D  angular 
extrapolation  algorithm  [11]  to  predict  the  high-frequency  RCS  of  complex  targets.  Our 
approach  is  based  on  modeling  the  induced  current  on  the  target  using  a  multiple-arrival 
model  that  closely  resembles  the  ray-optical  behavior  at  high  frequencies.  The  model 
coefficients  are  determined  by  the  superresolution  algorithm  ESPRIT  (Estimation  of 
Parameters  via  Rotation  Invariance  Technique)  [13.  14].  In  this  paper,  we  extend  our 
algorithm  to  two  dimensions  to  carry  out  simultaneous  frequency-aspect  extrapolation. 


We  begin  by  proposing  a  2-D  multiple-arrival  model  in  the  frequency-aspect  domain. 
Next,  we  implement  a  2-D  ESPRIT  algorithm  [15,  16]  to  estimate  the  model  coefficients 
from  a  limited  number  of  frequency-aspect  current  data  computed  using  a  fast  multipole 
solver.  Once  the  model  is  determined,  the  currents,  and  subsequently  the  RCS,  at  other 
frequencies  and  aspects  can  thus  be  computed. 

This  paper  is  organized  as  follows.  In  Section  II,  we  examine  the  frequency  and 
angular  dependency  of  the  current  phase  and  propose  a  multiple-arrival  model  in  the 
frequency-aspect  domain.  In  Section  III,  a  2-D  ESPRIT  algorithm  is  developed  and 
applied  to  estimate  the  time-of-arrival  and  cross  range  parameters  in  the  frequency-aspect 
model  from  a  few  calculated  frequency-aspect  samples.  In  Section  IV,  the  extrapolation 
algorithm  is  tested  on  a  canonical  2-D  cylinder-plate  structure.  The  result  is  compared 
against  both  the  exact  result  and  that  obtained  by  using  the  well-known  bistatic 
approximation  in  angle  [12,17,18].  The  extrapolated  results  are  also  generated  as  a 
function  of  the  number  of  calculated  samples  to  investigate  the  convergence  of  the 
algorithm.  Finally,  we  apply  the  algorithm  to  extrapolate  the  RCS  and  generate  a  2-D 
ISAR  image  of  the  VFY-218  airplane  at  UHF  band.  The  extrapolated  result  is  compared 
against  the  image  obtained  from  chamber  measurement  data  [19].  Conclusions  are  given 
in  Section  V. 

II.  2-D  MODEL  FOR  FREQUENCY-ASPECT  EXTRAPOLATION 

We  shall  first  formulate  a  physical  model  for  the  current  induced  on  the  surface  of 
a  target  due  to  an  incident  wave.  The  induced  current  at  any  point  is  in  general  the  result 
of  multiply  incident  waves  from  both  the  direct  excitation  and  multiple  scattering  from 


other  parts  of  the  target.  Therefore,  we  postulate  that  the  induced  current  can  be  written 
as  a  sum  of  multiply  incident  waves,  each  with  a  different  travel  paths,  as  shown  in  Fig. 
1.  If  we  denote  the  down  range  direction  with  respect  to  the  incident  wave  as  x  and  the 
cross  range  direction  as  y,  the  current  at  S  as  a  function  of  frequency / can  be  written  as 

K  -j=^-dk 

J(f,S)  =  ^Ake  c  ,dk  =  xk  +4  (1) 

k=  1 

where  K  is  the  number  of  incident  waves  arriving  at  S  and  c  is  the  speed  of  light  in  free 
space.  In  the  above  definition  of  the  path  length  dk,  we  let  (xk,  yk)  be  the  first  hit  point  on 
the  target  due  to  the  incident  wave,  and  lk  be  the  total  intermediate  path  length  of  the 
multiple  scattering  mechanism  from  the  first  hit  point  to  point  5.  Ak  is  the  amplitude 
coefficient  for  each  mechanism  and  is  assumed  to  be  frequency  independent.  Among  the 
three  mechanisms  illustrated  in  Fig.l,  mechanism  1  is  the  direct  incident  wave  from  the 
source.  Therefore,  l/=0  and  (xIt  yt)  corresponds  to  point  S.  Mechanisms  2  and  3  are 
respectively  a  once-scattered  and  a  twice-scattered  wave  before  arriving  at  5.  At  high 
frequencies,  this  model  is  expected  to  be  quite  sparse,  i.e.,  only  a  few  terms  are  needed  to 
adequately  describe  the  scattering  physics  at  an  arbitrary  point  S  on  the  target  surface. 
This  1-D  model  was  used  to  achieve  frequency  extrapolation  at  a  fixed  aspect  angle  in 
our  earlier  work  [10,12]. 

Angular  dependency  can  also  be  incorporated  in  the  above  model.  We  assume 
that  all  the  intermediate  scattering  points  and  the  amplitude  coefficient  for  each 
mechanism  remain  fixed  as  the  incident  angle  is  varied,  as  illustrated  in  Fig.l.  This 
assumption  was  found  to  be  fairly  accurate  for  ray-optical  fields  under  small  angular 
variation  [20].  We  have  also  applied  it  to  achieve  angular  extrapolation  at  a  fixed 


frequency  for  iterative  moment  solvers  [11].  When  we  combine  the  aspect  behavior  with 
the  frequency  model  in  (1),  we  arrive  at  the  following  model  for  the  current  as  a  function 
of  both  frequency  and  incident  angle: 


K 


J(f,d,S)  =  ^Ake 

k=\ 


-  j~~~~(xki  COS  &+>'kl  s'n  0+4i ) 
C 


(2) 


Equation  (2)  contains  three  unknowns  in  the  phase  function  of  each  mechanism.  Next  we 
use  the  small-angle  approximation  cos6~l  to  arrive  at  an  expression  with  two  unknowns 
in  the  phase: 


J{f,d,S)  =  ^Ake  ‘ 


K  (</t,cos0+vt,.sin(?)  ^  ,  -j(k,du+k,yti) 


=X-v 


(3) 


k= 1 


k= 1 


where kx  =  2nfcos6/c  and  k  =  2nfsin6/c.  Note  that  equation  (3)  can  be  further 
approximated  to  completely  decouple  the  frequency  and  aspect  variable: 


J(f,9,S)  =  ^Ake 
k= 1 


r  .  2nfc  .  2xf 

K  -j — —yki0  ~J—dki 


(4) 


if  we  use  sin0~6  and  replace  the  frequency  variable  in  the  first  exponential  by  the  center 
frequency  fc.  Equation  (4)  then  reduces  to  equation  (3b)  in  [11]  when  only  aspect 
variation  is  considered.  For  the  2-D  extrapolation  in  this  work,  we  choose  to  use  the 
model  in  (3)  since  it  is  more  accurate  when  the  aspect  range  is  large.  Equation  (3)  is  in 
the  form  of  a  sum-of-exponential  model,  with  linear  phase  dependence  with  respect  to 
both  kx  and  fcv„  Consequently,  the  2-D  superresolution  algorithm  ESPRIT  [17]  can  be 


applied  to  equally  spaced  kx  and  ky  data  samples.  Next,  we  shall  describe  the 
implementation  of  the  2-D  ESPRIT  algorithm  to  estimate  the  unknown  parameters  in  the 


model. 


III.  TWO-DIMENSIONAL  ESPRIT  ALGORITHM 


In  [16],  a  2-D  ESPRIT  algorithm  was  developed  for  estimating  the  direction-of- 
arrival  in  2-D  antenna  array  problems.  Here,  we  shall  show  that  with  some  minor 
modifications,  it  can  be  applied  to  estimate  the  parameters  A* ,  dk  and  yk  in  equation  (3). 
First,  we  assume  that  the  parameters  are  to  be  estimated  from  known  current  values 
solved  at  M  x  N  equally  spaced  samples  in  the  kx  -  k  plane,  where  M  is  the  number  of 


samples  in  kx  and  N  is  the  number  of  samples  in  k v.  These  M  x  N  samples  are  similar  to 
the  elements  of  a  2-D  antenna  array  described  in  [16].  We  define: 


pk=e-^‘, 


(5) 


where  Akx and  Aky  are  the  sampling  intervals  inland  k,  respectively.  If  we  shift  the 
origin  of  the  variables  kx  and  k  to  zero,  we  can  rewrite  equation  (3)  into  a  form  similar 
to  equation  (9)  of  [16]: 

z«  =i^pr'<?r  +na,  i  =  1,2. . M.  /  =  1.2 . N  (6) 

*=1 

where  zu  is  the  current  at  the  /th  kx  value  and  the  /th  k.  value.  sk  is  the  modified 

amplitude  coefficient  for  the  Ath  mechanism  to  account  for  the  origin  shift,  n is  assumed 

to  be  white  Gaussian  noise,  which  in  our  case  is  used  to  model  the  numerical  error  in  the 
current  computation.  Since  ESPRIT  postulates  such  a  sum-of-exponential  model  with 
additive  white  noise,  an  averaging  procedure  has  to  be  performed  to  smooth  out  the  noise 
to  obtain  the  correct  estimate  of  the  covariance  matrix.  In  [16],  the  averaging  is 
performed  naturally  in  the  time  domain.  However,  for  the  problem  at  hand,  the  time 


dimension  does  not  exist  and  has  to  be  synthesized.  This  can  be  accomplished  via  a  sub- 
array  processing  technique.  As  shown  in  Fig.  2,  a  sub-array  size  is  chosen  to  be  half  of 
the  original  array.  Then  shifting  the  sub-array  one  data  sample  at  a  time  in  either  kx  or 

ky  will  result  in  a  new  sub-array,  which  can  be  considered  as  the  array  at  a  new  time 
index.  If  we  assume  M  and  N  are  even  numbers,  the  total  number  of  sub-arrays  that  can 


be  generated  in  this  manner  is(y-  +  l)(y  + 1).  The  signal  vectors  z(t),  s(t)  and  n(t) 


are  defined  as  the  reshaped  column  vectors  of  the  rth  sub-array  as  follows: 


z(t)  =  U,  1  (0-  •  •  zlL  (0  ■  ■  ■  zn  (0- •  *  ZIL  it  )]r 

(7) 

s(t)  =  [sl(t)  s2(t)---sK(t))T 

(8) 

(9) 

r  M  J  r  N 
In  the  above  expressions, I  =~  and  L  =  — , 


which  are  the  maximum  row  and  column 


indices  of  each  sub-array.  The  array  covariance  matrix  can  thus  be  defined  exactly  like 
formula  (15)  in  [16].  Following  the  procedures  in  [16],  we  estimate  the  parameters  pk 
and  qk.  The  time-of-arrival  parameter  dk  and  the  cross-range  parameter  yk  are  obtained 
using: 


M, 


-^Pk’ 


>'k  = 


M, 


-^Qk 


(10) 


Finally,  the  amplitude  coefficient  Ak  can  be  generated  by  solving  equation  (3)  via  a  least 
squares  procedure.  Noted  that  the  maximum  order  number  K  permitted  by  the  algorithm 
is: 


£<min[(4-l)4,(4-l)4]-l 


(ID 


For  the  scattering  problems  we  have  examined,  K  is  usually  chosen  to  be  3  or  4.  The 
minimum  number  of  frequency  and  angular  samples  can  then  be  selected  to  satisfy  the 
above  criterion. 

After  the  coefficients  Ak ,  dk  and  y*  are  estimated,  we  assume  equation  (3)  also 
holds  for  other  frequencies  and  aspects.  Therefore,  the  induced  current  can  be 
extrapolated  to  other  kx  and  ky  values  of  interest.  An  IS  AR  image  is  then  generated  using 

standard  Fourier  processing  of  the  kx-k>,  data.  To  summarize,  the  extrapolation 

procedure  is  carried  out  by  first  selecting  a  number  of  densely  sampled  points  in  a  limited 
frequency-aspect  range  and  solving  the  scattering  problems  at  these  points.  Usually  the 
computed  frequencies  are  chosen  close  to  the  low  frequency  end  and  the  aspect  angles  are 
centered  about  the  central  angle  of  interest.  Second,  the  2-D  ESPRIT  algorithm  is  applied 
to  estimate  the  model  parameters  of  the  current  at  each  point  on  the  target  surface.  Third, 
the  induced  current  is  extrapolated  to  a  wider  kx  and  k-  range  based  on  the  model 

coefficients  generated  from  2-D  ESPRIT  processing.  Finally,  the  current  is  integrated  to 
generate  the  scattered  far  field  as  a  function  of  frequency  and  aspect  and  the  results  are 
used  to  generate  the  desired  ISAR  image. 

IV.  NUMERICAL  RESULTS 

To  validate  the  frequency-aspect  extrapolation  algorithm,  we  first  consider  a  2-D 
cylinder-plate  structure  shown  in  Fig.  3.  The  diameter  of  the  cylinder  is  4.2  m  and  the 
length  of  the  plate  is  20  m.  The  origin  of  the  cylinder  and  the  center  of  the  plate  are 
separated  by  6.2  m.  The  frequency  band  of  interest  is  from  0.3  GHz  to  0.65  GHz  and  the 
observation  angle  is  from  25°  to  65°.  Fig.  4  shows  the  reference  ISAR  image  generated 


from  71x  81=5751  computed  points  in  the  frequency-aspect  plane.  The  image  has  a 
dynamic  range  of  40  dB.  In  this  structure,  strong  multiple  scattering  mechanisms 
between  the  cylinder  and  the  plate  dominate  the  backscattering.  We  can  see  in  the  image 
the  features  corresponding  to  the  direct  scattering  from  the  cylinder  (labeled  as  (i)),  the 
front  edge  point  of  the  plate  (ii),  and  the  shadow  boundary  cast  on  the  plate  by  the 
cylinder  (iii).  Additionally,  there  are  other  range-delayed  features  corresponding  to  the 
multiple  scattering  mechanisms.  For  comparison,  we  calculate  the  current  at  71 
frequency  samples  for  one  aspect  angle,  and  use  the  well-known  bistatic  approximation 
[17,  18]  to  extrapolate  the  RCS  to  other  aspect  angles.  Fig.  5  plots  the  resulting  ISAR 
image  from  the  extrapolated  frequency-aspect  data  using  the  bistatic  approximation.  In 
this  image,  only  the  direct  scattering  features  ((i)  and  (ii))  are  correctly  predicted  while 
the  higher-order  scattering  features  are  poorly  predicted  in  either  position  or  amplitude. 
This  is  because  the  bistatic  approximation  is  based  on  physical  optics,  in  which  the 
current  is  assumed  to  be  excited  by  only  the  direct  incident  wave. 

We  next  construct  the  ISAR  image  using  the  2-D  frequency-aspect  extrapolation 
algorithm.  We  choose  8x8  points  at  the  low  frequency  end  and  about  the  central  angle 
of  the  original  7 lx  81  points.  After  the  time-of-arrival  and  cross-range  parameters  are 
extracted  by  using  2-D  ESPRIT,  the  induced  currents  and  scattered  far  fields  are 
extrapolated  to  the  kx-kx  aperture  of  the  original  samples.  In  this  manner,  the  ISAR 

image  can  be  generated  directly  by  a  2-D  FFT  of  the  kx-k}.  data  without  the  polar 

reformatting  operation  [2],  The  resulting  ISAR  image  is  plotted  in  Fig.  6(a).  As  we  can 
see,  most  of  the  features  are  predicted  in  the  correct  position.  However,  some  of  the 
weaker  features  are  highly  defocused.  We  repeat  this  process  using  2-D  ESPRIT 


extrapolation  for  9  x  9,  10  x  10  and  11x11  input  points.  The  resulting  images  are 
shown  in  Figs.  6(b),  6(c)  and  6(d),  respectively.  From  this  series  of  images,  we  see  that 
the  strong  scattering  features  are  quite  stable  throughout  while  the  weaker  features  begin 
to  converge  as  the  number  of  input  points  is  increased.  We  find  that  the  correlation 
indices  between  the  extrapolated  images  and  the  reference  image  are  respectively  0.72, 
0.80,  0.82  and  0.88.  This  shows  a  steady  improvement  as  the  number  of  input  points  is 
increased.  For  the  case  of  11  x  11  input  points,  nearly  all  of  the  features  in  the 
extrapolated  image  agree  well  with  those  in  the  reference  image  from  the  brute-force 
calculation.  This  corresponds  to  a  7:1  extrapolation  ratio  in  each  dimension. 

Next,  we  use  the  2-D  extrapolation  algorithm  to  simulate  the  ISAR  image  of  the 
benchmark  VFY-218  airplane.  The  simulation  is  carried  out  using  the  multi-level  fast 
multipole  code  FISC  [21]  on  a  Pentium  II  450  MHz  PC.  The  total  computation  time 
would  take  about  200  days  if  a  brute-force  calculation  is  carried  out  over  the  required 
41x41  frequency-aspect  samples  (to  achieve  0.5  meter  resolution  in  range  and  cross 
range).  Based  on  the  7: 1  criterion  found  above,  the  actual  computation  is  only  carried  out 
for  6  x  6  =  36  points  from  267  MHz  to  297  MHz  in  frequency  and  127°  to  132°  in  aspect. 
The  total  computation  time  is  40  hours.  Note  that  the  time  savings  over  the  brute-force 
calculation  is  even  greater  than  the  (41/6)2  ratio.  This  is  because  these  36  points  are 
chosen  at  the  low  end  of  the  frequency  band,  and  take  far  less  time  to  compute  than  those 
at  the  high  frequency  end.  The  induced  current  and  far  fields  are  then  extrapolated  to  a 
kx-k>.  grid  from  267  MHz  to  533  MHz  in  frequency  and  from  1 10°  to  140°  in  aspect.  The 

resulting  ISAR  image  of  the  airplane  at  130°  is  shown  in  Fig.  7.  Fig.  8  shows  the  ISAR 
image  constructed  from  chamber  measurement  data  at  the  same  look  angle  [19]. 


Comparing  these  two  images,  we  observe  that  nearly  all  of  the  key  features  in  the 
measurement  image  are  predicted  in  the  simulated  image  from  using  FISC  and  the  2-D 
extrapolation  algorithm.  As  expected,  the  image  in  Fig.  7  is  much  cleaner  than  the 
previously  predicted  image  shown  in  Fig.  11  of  [12],  which  was  simulated  from  the  1-D 
frequency  extrapolation  procedure  in  combination  with  the  bistatic  approximation.  The 
improvement  is  due  to  the  fact  that  the  multiple  scattering  mechanisms  are  correctly 
incorporated  in  both  dimensions  in  our  new  model.  There  still  exist  some  low  dynamic 
range  noises  in  the  predicted  image,  which  can  be  improved  with  more  computed  points. 

V.  CONCLUSIONS 

A  2-D  frequency-aspect  extrapolation  algorithm  has  been  developed  to  accelerate 
ISAR  image  simulation  using  fast  multipole  solvers.  A  2-D  multiple-arrival  model 
consistent  with  high-frequency  scattering  physics  was  proposed  to  model  the  induced 
currents  on  the  target.  A  2-D  ESPRIT  algorithm  was  developed  to  estimate  the  model 
parameters  from  a  limited  number  of  computed  data  samples  in  frequency  and  aspect. 
The  model  was  then  extrapolated  to  other  frequencies  and  aspects  to  arrive  at  broadband, 
wide-angle  RCS  data  for  ISAR  image  construction.  This  algorithm  has  been  tested  using 
a  canonical  cylinder-plate  structure  to  evaluate  its  performance.  It  was  found  that  a  7:1 
extrapolation  ratio  in  each  dimension  can  been  achieved.  Furthermore,  the  input  data 
samples  can  be  chosen  in  the  low  frequency  range  for  greater  computational  payoff. 
Lastly,  the  ISAR  image  of  the  benchmark  VFY-218  airplane  at  UHF  band  has  been 
predicted  using  the  fast  multipole  solver  FISC  and  the  2-D  extrapolation  algorithm.  The 
resulting  image  compared  favorably  with  that  obtained  from  chamber  measurement  data. 


The  computation  time  savings  over  the  brute-force  calculation  was  about  two  orders  of 
magnitude. 
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Fig.  1 .  Multiple-arrival  model  for  the  induced  current. 
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Fig.  4.  Reference  ISAR  image  generated  using  the  brute-force 
MoM  calculations  at  71x81=5751  points. 


<D  30 

to 

c 

ca 


CO 

CO 


25 


O 

U  20 


“s 

m- 

fc 

i 

a 

125 


-5 


20  25 

range 


(dBm) 


Fig.  5.  ISAR  image  generated  from  71  calculated  frequency 
points  and  the  bistatic  approximation. 
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Fig.  7.  ISAR  image  of  the  VFY-218  at  130  degrees  from  nose-on 
generated  from  the  extrapolated  result  using  2-D  ESPRIT 
with  36  FISC-computed  points. 
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Fig.  8.  ISAR  image  of  the  VFY-218  at  130  degrees  from  nose-on 
generated  from  chamber  measurement  data  [19]. 
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Abstract 

A  frequency  extrapolation  technique  is  proposed  to  obtain  the  broad  band  radiation 
patterns  of  antenna-platform  radiation  problems.  A  frequency-dependent  time-of-arrival 
model  is  applied  to  the  induced  current  computed  at  low  frequencies.  The  model 
parameters  are  estimated  by  applying  the  ESPRIT  superresolution  algorithm  to  the 
computed  data.  A  pre-multiplication  scheme  in  conjunction  with  the  complex  time-of- 
arrival  estimation  from  ESPRIT  is  used  to  determine  the  additional  frequency-dependent 
factors  in  the  model.  The  current  at  high  frequencies  is  then  extrapolated  based  on  the 
model  and  the  radiated  field  is  computed  using  the  extrapolated  current.  The  algorithm  is 
applied  to  several  2D  and  3D  antenna-platform  radiation  problems.  The  extrapolated 
radiation  patterns  are  found  to  be  in  good  agreement  with  those  generated  by  brute  force 
computation.  Since  the  current  required  for  the  extrapolation  is  only  computed  at  lower 
frequencies,  large  savings  in  computation  time  and  memory  can  be  achieved. 


I.  Introduction 


The  numerical  solution  of  electromagnetic  scattering  and  radiation  problems  from 
complex  structures  is  usually  very  computer  intensive  in  time  and  memory,  especially  at 
high  frequencies.  In  the  standard  method  of  moments  (MoM),  for  example,  the 
computation  complexity  scales  as/6  where  /is  the  frequency.  Even  for  fast  solvers  such 
as  the  fast  multipole  method  [1,2],  the  frequency  scaling  is  still  between  f2  and /3.  Thus 
it  is  desirable  to  obtain  the  broad  band  response  from  the  computation  results  at  only  a 
few  frequencies.  This  can  be  accomplished  by  fitting  the  computed  frequency  response  to 
a  reduced-order  model  and  extracting  the  model  parameters  from  the  data.  Approaches 
based  on  model-based  parameter  estimation  have  been  studied  extensively  and  applied  to 
many  aspects  in  electromagnetic  computation  [3-11].  The  appropriate  reduced-order 
models  to  parameterize  the  physical  observables  are  different  depending  on  the  frequency 
region  of  interest.  In  the  low  frequency  range,  the  rational  function  model  is  widely  used, 
while  in  the  high  frequency  range,  the  time-of-arrival  model  is  the  preferred  choice  in 
characterizing  the  ray-optical  behaviors  of  fields  and  currents. 

In  this  paper,  we  address  the  model-based  frequency  extrapolation  of  antenna- 
platform  radiation  problems.  It  is  well  known  that  the  radiation  patterns  of  antennas  are 
strongly  affected  by  the  mounting  platform.  Since  the  simulation  of  such  radiation 
problems  is  very  costly  to  generate  as  a  function  of  frequency,  we  set  out  to  extrapolate 
the  high  frequency  response  from  computed  data  at  lower  frequencies.  Our  approach 
entails  fitting  the  computed  current  at  low  frequencies  to  a  time-of-arrival  model  and 
estimating  the  model  parameters  by  the  superresolution  algorithm  ESPRIT  [12]. 
Previously,  this  approach  has  been  applied  to  radar  signature  prediction  with  good 


success  [10].  In  that  work,  the  coefficients  of  the  time-of-arrival  model  were  assumed  to 
be  frequency  independent.  In  the  radiation  problem,  the  primary  radiation  from  the 
antenna  is  usually  frequency  dependent.  Furthermore,  when  there  exist  higher-order 
multiple  interactions  on  the  platform,  the  amplitude  of  each  mechanism  is  in  general  also 
frequency  dependent  [13-17].  Thus  a  more  accurate  model  is  needed  to  parameterize  the 
current.  Here  we  generalize  the  extrapolation  algorithm  by  using  a  frequency-dependent 
time-of-arrival  model.  To  extract  the  additional  frequency  dependence  in  this  model,  we 
adopt  an  approach  similar  to  the  one  proposed  in  [15],  which  was  developed  to 
effectively  extract  the  frequency  dependency  of  the  scattering  mechanisms  in  measured 
backscattered  data.  Our  results  show  that  the  radiation  pattern  can  be  accurately 
extrapolated  based  on  the  frequency  dependent  model. 

This  paper  is  organized  as  follows.  Section  II  describes  the  frequency-dependent 
time-of-arrival  model  and  the  procedure  to  determine  the  model  parameters.  Section  III 
gives  a  numerical  example  of  the  extrapolation  algorithm  using  simulated  data.  The 
performance  of  the  algorithm  in  the  presence  of  noise  is  investigated  and  the  errors  in  the 
estimation  of  model  parameters  are  quantified.  In  Section  IV  we  apply  the  algorithm  to 
extrapolate  the  induced  current  in  antenna-platform  radiation  problems.  The  radiation 
patterns  predicted  at  high  frequencies  are  compared  against  the  reference  MoM 
computation  in  both  2D  and  3D  platforms.  Conclusions  are  given  in  Section  V. 

II.  Frequency  Dependent  Model  and  Determination  of  Model  Parameters 

We  postulate  that  the  current  induced  on  a  complex  platform  due  to  illumination  from 
a  primary  source  can  be  well  described  by  a  time-of-arrival  model  at  high  frequencies. 


The  different  time-of-arrival  terms  correspond  to  the  different  incident  wave  mechanisms 
from  both  the  direct  antenna  radiation  and  higher-order  scattering  from  other  parts  of  the 
platform,  as  illustrated  in  Fig.  1.  Therefore,  the  current  at  an  arbitrary  point  P  on  the 
platform  can  be  expressed  as: 

Jm  =  ^ane-j^  (1) 

n 

where  tn  is  the  arrival  time  of  the  nth  incident  wave  and  an  is  its  amplitude.  In  [10],  an 

was  assumed  to  be  independent  of  frequency.  However,  in  the  antenna-platform  radiation 
problem,  the  primary  radiation  from  the  antenna  is  in  general  frequency  dependent.  For 
instance,  the  radiation  strength  of  a  dipole  antenna  is  proportional  to  the  co  (or  co'12  for  a 
line  source  in  2D  problems).  Furthermore,  for  complex  platforms,  the  incident  waves  to  a 
specific  current  element  can  come  from  the  scattered  fields  from  other  parts  of  the 
platform.  These  secondary  incident  fields  are  also  frequency  dependent.  For  canonical 
shapes,  their  exact  frequency  dependencies  are  predicted  by  the  geometrical  theory  of 
diffraction  (GTD),  and  are  in  the  form  of  of  where  a  is  the  frequency  factor.  For 
example,  the  strength  of  the  scattered  field  from  comer  diffraction  has  a  frequency 
dependency  of  co1,  and  that  from  edge  diffraction  has  a  frequency  dependency  of  coU2. 
Therefore,  the  time-of-arrival  model  can  be  improved  by  incorporating  the  frequency 
factor: 

(2) 

n 

where  an  is  the  frequency  factor  for  each  incident  wave.  Compared  with  (1),  this  model  is 
better  matched  to  the  actual  physics  of  the  radiation  and  scattering  mechanisms. 


Since  we  are  attempting  to  extrapolate  the  frequency  response  to  a  much  broader 
range,  the  accurate  estimation  of  the  frequency  factors  is  critical.  A  small  error  in  a  will 
result  in  dramatic  difference  in  amplitude  at  frequencies  in  the  extrapolated  region. 
However,  the  existing  superresolution  algorithms  based  on  eigenspace  decomposition 
(e.g.  ESPRIT  and  MUSIC)  apply  only  to  signals  in  the  form  of  (1).  For  instance,  the 
ESPRIT  algorithm  [12]  estimates  ( an ,  t„)  from  uniformly  sampled  current  data  from  a>\  to 
com  based  on  the  data  model: 

)  =  £  «* «  +  n«°m )’  <ym  =  CO, ,  0)2 , ' •  •  ■ • ,  coM  (3) 

n 

where  n  denotes  additive  Gaussian  white  noise.  This  model  can  be  easily  extended  to 
estimate  complex  time-of-arrivals  via  analytic  continuation.  If  we  separate  the  real  and 
imaginary  part  of  f,„  the  resulting  model  can  be  written  as 

J(co)  =  ^anew  lm<r'  V'"  Re(F"  >  (4) 

n 

where  the  tilda  symbol  is  used  to  indicate  a  complex  number.  This  is  the  well-known 
complex  exponential  model  and  has  been  used  to  approximately  model  the  frequency 
dependence  in  an  [13].  Comparing  (4)  to  the  frequency-dependent  time-of-arrival  model 
(2),  we  find  that  the  only  difference  between  the  two  is  the  frequency  dependency,  which 

is  in  the  form  of  of  in  (2)  and  eQ)lm('1')  in  (4).  For  a  narrow  band  signal,  (4)  is  usually  a 
good  approximation  of  (2),  and  can  be  simply  used  to  model  the  frequency  dependence 
within  this  band.  However,  in  order  to  achieve  accurate  extrapolation  results,  the  complex 
exponential  model  is  a  poor  choice  since  it  is  not  properly  matched  to  the  underlying 
physics. 


To  better  determine  the  frequency  factor  a„,  we  adopt  an  approach  similar  to  the  one 
proposed  in  [15]  based  on  a  pre-multiplication  concept.  We  multiply  the  data  by  an 
assumed  frequency  factor  (0~a” .  The  data  model  then  becomes 

Jfco)  =  c )  =  ^antyla”'Q'”)e"-'""  (5) 

n 

Next,  we  apply  ESPRIT  and  the  complex  exponential  model  on  J{co)  to  estimate  the  real 
and  imaginary  parts  of  tn .  Note  that  if  the  frequency  dependence  of  the  pre-multiplier  is 

exactly  matched  to  that  of  a  particular  mechanism,  the  resulting  exponential  term  will 
have  no  frequency  dependence.  Consequently,  the  estimated  complex  time-of-arrival 
term  will  have  a  zero  imaginary  part.  Therefore,  by  repeatedly  pre-multiplying  the  signal 
using  different  values  of  am  and  observing  the  imaginary  parts  of  the  resulting  7n ,  we  can 

detect  the  correct  frequency  dependence  an  whenever  Imf  7n )  goes  to  zero.  The  implicit 
assumption  of  this  approach  is  that  the  mismatched  terms  will  not  significantly  distort  the 
estimation  of  Imf  7n)  of  the  matched  term.  This  is  usually  true  when  the  arrival  times 

Ref  7n)  are  well  separated.  When  two  or  more  of  the  arrival  times  are  very  close  to  each 
other,  this  algorithm  becomes  less  accurate  due  to  the  interference  between  the 
exponential  terms.  The  imaginary  part  of  7n  may  not  vanish  even  when  the  corresponding 

frequency  factor  is  matched  by  the  pre-multiplication.  Quantitative  evaluation  of  the 
resulting  error  will  be  discussed  in  detail  in  the  next  section.  Once  an  and  tn  are  estimated, 
the  amplitude  cin  can  be  easily  found  by  the  standard  least  squares  procedure. 

III.  Algorithm  Testing  and  Error  Estimation  Using  Simulated  Data 

To  test  the  extrapolation  algorithm,  we  consider  the  signal 


J(0))  =  (1  +  0.4  j)(t)-U2e-J2S<a  +  (4.3  -  2.2  yV' V'1 3<a  +  n,  co  =  5,6,...,35  (6) 

where  n  is  additive  Gaussian  white  noise.  The  signal-to-noise  ratio  (SNR)  is  set  to  10  dB. 
We  now  use  a  two-term  time-of-arrival  model  to  extrapolate  the  signal  based  on  its  first 
10  samples.  The  actual  signal  is  plotted  as  the  solid  line  in  Fig.  2.  The  extrapolation 
results  of  three  different  models  are  plotted.  The  frequency  independent  model,  which  is 
plotted  in  dashed  dot,  is  obtained  by  the  model  in  (1)  with  real  time-of-arrival  tn 
determined  by  ESPRIT.  This  model  matches  badly  with  the  original  curve,  even  in  the 
sampled  region.  The  complex  exponential  model,  plotted  in  dashed  line,  uses  the  model 
in  (4)  with  complex  tn .  The  resulting  curve  is  in  good  agreement  with  the  actual  signal 

within  the  sampled  region,  but  not  in  the  extrapolated  region.  This  indicates  that  although 
the  model  (4)  can  be  a  very  good  approximation  to  the  actual  signal  in  the  sampled 
region,  it  cannot  be  used  to  accurately  extrapolate  the  signal  because  it  does  not  have  the 
correct  frequency  dependency.  The  frequency-dependent  time-of-arrival  model,  which  is 
plotted  in  dots,  matches  the  best  with  the  original  signal  in  both  the  sampled  region  and 
the  extrapolated  region.  Its  parameters  are  estimated  using  the  method  described  in  Sec. 
II. 

To  show  the  detail  workings  of  the  pre-multiplier  procedure,  the  imaginary  parts  of  Tn 

corresponding  to  the  two  terms  in  (6)  are  plotted  as  functions  of  am  in  Figs.  3(a)  and  3(b) 
at  the  SNR  levels  of  and  10  dB,  respectively.  The  two  curves  cross  zero  at  the 
corresponding  frequency  factor  of  -0.5  and  0.5.  It  can  be  observed  that  the  shapes  of  the 
curves  change  when  the  SNR  is  decreased,  especially  for  the  term  with  smaller  power 
(a^-O.S).  However,  the  positions  of  the  zero-crossing  points  vary  only  weakly  as  the 


SNR  is  decreased. 


Next,  we  examine  the  error  performance  of  the  extrapolation  algorithm  through 
numerical  simulation.  Among  the  three  sets  of  model  parameters  an,  and  tn,  the  time- 
of-arrival  tn  is  directly  estimated  by  the  ESPRIT  algorithm.  When  frequency  dependent 
mechanisms  are  present  in  the  signal,  the  estimation  error  in  tn  is  expected  to  be  larger 
due  to  the  interference  between  the  exponential  terms,  especially  when  two  or  more 
arrival  times  are  very  close  to  each  other.  In  this  example,  we  estimate  the  arrival  times  in 
a  simulated  signal  consisting  of  two  frequency  dependent  exponential  terms  using 
ESPRIT.  As  a  comparison,  we  also  determine  the  arrival  times  in  a  two-term  frequency 
independent  signal.  In  both  signals,  the  power  of  the  stronger  term  is  five  times  that  of  the 
weaker  term.  For  the  frequency  dependent  signal,  the  two  frequency  factors  are  set  to  0.5 
and  -0.5  for  the  stronger  and  weaker  terms,  respectively.  The  test  is  performed  1000 
times  subject  to  random  noise.  In  Fig.  4(a),  the  root  mean  squared  (RMS)  error  on  the 
time-of-arrival  estimate  of  the  weaker  term  is  plotted  as  a  function  of  the  separation 
between  the  two  arrival  times  and  for  different  SNR  conditions.  The  dashed  lines  are  the 
estimation  errors  from  the  frequency  independent  signal,  while  the  solid  lines  are  the 
errors  from  the  frequency  dependent  signal.  Both  axes  are  calibrated  in  terms  of  the 
Fourier  resolution,  which  is  the  reciprocal  of  the  available  data  bandwidth.  It  is  shown 
that  the  arrival  time  error  is  larger  for  the  frequency  dependent  signal,  especially  when 
the  SNR  is  high.  However,  the  degradation  is  not  too  severe  overall.  Good  estimation 
(with  error  less  than  0. 1  Fourier  bin)  is  still  achievable  for  signal  separation  well  within 
one  Fourier  bin. 

The  RMS  error  of  estimating  the  frequency  factor  an  of  the  weaker  signal  is  plotted  in 
Fig.  4(b)  for  the  frequency  dependent  signals  used  in  the  above  test.  We  observe  that  the 


error  decreases  as  the  separation  between  the  exponential  terms  increases  and  as  the  SNR 
increases,  similar  to  the  trends  in  Fig.  4(a).  When  the  two  arrival  times  are  close,  the 
interaction  between  the  two  terms  results  in  error  on  the  estimation  of  a„  and  the  error  is 
strongly  affected  by  the  noise  level.  We  can  further  decrease  the  error  by  imposing  the 
constraint  that  the  frequency  factors  are  integer  multiples  of  1/2,  as  dictated  by  GTD. 
This  is  implemented  when  we  deal  with  actual  electromagnetic  data  in  the  next  section. 

Finally  we  look  at  the  extrapolation  error  of  the  entire  frequency-dependent  time-of- 
arrival  model.  An  extrapolation  ratio  of  3  to  1  is  used  in  the  example,  i.e.  the  frequency 
response  is  extrapolated  to  three  times  the  original  bandwidth.  The  percentage  RMS  error 
is  plotted  in  Fig.  4(c).  When  the  separation  between  the  two  terms  is  large,  the 
extrapolation  error  behaves  similar  to  the  errors  of  tn  and  an.  When  the  two  terms  become 
too  close  in  time,  we  note  that  the  error  curves  actually  reach  a  maximum  and  then  drop 
off.  This  is  because  below  this  point,  the  two  exponential  terms  are  too  close  to  be 
resolved  within  the  extrapolation  region  of  interest.  For  larger  extrapolation  bandwidth, 
the  position  of  the  peak  should  move  toward  zero. 

IV.  Numerical  Results 

We  now  apply  the  frequency  extrapolation  technique  to  computation  data  from 
antenna-platform  radiation  problems.  In  the  first  example,  a  two-dimensional  structure 
shown  in  Fig.  5  is  considered.  The  platform  is  14m  in  length  and  3m  in  height.  The 
antenna  is  a  horizontal  line  source  placed  at  5m  above  the  platform.  The  induced  current 
on  the  platform  is  computed  from  0.1  to  0.5  GHz  at  21  frequency  points  using  2D  MoM. 
The  current  is  extrapolated  to  1.3  GHz  and  radiated  field  is  then  computed  based  on  the 


extrapolated  current.  Both  the  frequency-independent  and  frequency-dependent  time-of- 
arrival  models  are  used  to  perform  the  extrapolation,  and  the  resulting  radiated  fields  at 
the  observation  angle  of  6  =  40°  are  plotted  in  Figs.  6(a)  and  6(b),  respectively.  Also 
plotted  is  the  reference  MoM  results  obtained  via  brute  force  computation.  The  primary 
radiation  of  the  dipole  antenna  is  not  included  in  the  plots  so  that  we  can  better  observe 
the  secondary  radiation  from  the  platform.  It  is  obvious  that  the  frequency  dependency  in 
the  field  response  is  not  captured  by  the  frequency-independent  time-of-arrival  model, 
while  the  field  predicted  by  the  frequency  dependent  model  is  in  good  agreement  with  the 
computed  result.  The  time  domain  response  corresponding  to  Fig.  6(b)  is  plotted  in  Fig. 
6(c).  It  is  shown  that  the  time-domain  peaks  are  well  characterized  by  the  extrapolated 
field.  Finally,  the  radiation  patterns  obtained  through  brute  force  MoM  computation  and 
frequency  extrapolation  are  plotted  as  functions  of  frequency  and  observation  angles  in 
Figs.  7(a)  and  7(b),  respectively.  Very  good  agreement  is  observed  between  the  two 
figures.  For  a  quantitative  comparison,  the  correlation  index  between  the  two  figures  is 
computed  using  the  definition 

ff  E;{f.d)E,{f,d)dfdd 

R  =  - — u - ; -  (7) 

,(f,ef+\E2(f,ef)dfde 

where  E /  and  £?  are  the  radiated  fields  in  linear  scale  obtained  by  extrapolation  and  MoM 
computation,  respectively.  The  correlation  index  between  the  two  figures  is  0.9992  in  the 
sampled  region  and  0.9857  in  the  extrapolated  region.  As  expected,  the  correlation  is 
lower  in  the  extrapolated  region.  However,  the  drop  off  is  relatively  small. 

Next,  we  look  at  a  3D  platform  shown  in  Fig.  8.  The  antenna  is  a  horizontal  dipole 
oriented  in  the  x  direction.  The  induced  current  is  computed  from  0.1  to  0.36  GHz  at  13 


frequencies  and  extrapolated  to  0.7GHz.  The  computation  is  carried  out  using  FISC  [18], 
which  is  a  3D  MoM  code  based  on  the  fast  multipole  method.  The  extrapolated  frequency 
and  time  domain  responses  at  the  observation  angle  <f>et  =  30°,  <paz  =  -60°  are  plotted  in 
Figs.  9(a)  and  9(b),  respectively.  Also  plotted  for  comparison  are  the  reference  responses 
computed  by  FISC  via  brute  force.  The  major  radiation  features  are  captured  in  both 
domains  by  the  extrapolation.  In  the  time  domain  response,  we  observe  that  the 
amplitudes  of  some  of  the  weaker  peaks  are  not  correctly  predicted  by  the  extrapolation, 
although  their  arrival  times  are  estimated  quite  accurately.  We  believe  this  is  due  to  the 
estimation  error  of  Or,  in  the  current  parameterization,  as  we  expect  larger  errors  in  the 
frequency  factors  for  the  weaker  time-of-arrival  terms.  Figs.  10(a)  and  10(b)  show  the 
reference  and  extrapolated  radiation  patterns  as  functions  of  frequency  and  azimuth  angle 
when  the  elevation  angle  is  fixed  at  50°.  Good  qualitative  agreement  is  observed.  The 
correlation  index  between  the  two  figures  is  found  to  be  0.9980  in  the  sampled  region  and 
0.9742  in  the  extrapolated  region,  respectively.  This  result  is  a  little  worse  than  the  2D 
example  because  the  noise  level  of  the  FISC-computed  current  is  higher  than  that  of  the 
2D  MoM  code.  As  was  shown  in  the  last  section,  the  extrapolation  performance  is 
affected  by  the  SNR  in  the  computed  data.  The  computation  time  of  the  brute  force 
reference  results  is  about  50  hours  on  a  Pentiumll  400MHz  PC,  while  the  total 
computation  time  to  carry  out  the  EM  computation  in  the  sampled  region  and  the 
extrapolation  process  is  7  hours. 


V.  Conclusions 


In  this  paper,  the  frequency-dependent  time-of-arrival  model  has  been  applied  to 
extrapolate  the  induced  current  and  radiation  pattern  in  antenna-platform  radiation 
problems.  The  model  parameters  are  estimated  by  applying  the  ESPRIT  superresolution 
algorithm  to  the  computed  data.  A  pre-multiplication  scheme  in  conjunction  with  the 
complex  time-of-arrival  estimation  from  ESPRIT  is  used  to  determine  the  additional 
frequency-dependent  factors.  The  performance  of  the  algorithm  in  the  presence  of  noise 
has  been  evaluated  based  on  simulated  data  and  errors  in  the  estimation  of  model 
parameters  have  been  quantified.  Our  results  show  that  the  method  is  quite  robust.  The 
algorithm  has  been  applied  to  extrapolate  the  induced  currents  and  radiation  patterns  in 
both  2D  and  3D  antenna-platform  radiation  problems.  The  radiation  patterns  computed 
from  the  extrapolated  currents  have  been  found  to  be  in  good  agreement  with  those 
generated  by  brute  force  computation. 

Although  the  determination  of  model  parameters  is  more  complicated,  the  frequency 
dependent  model  show  significant  performance  improvement  over  the  frequency 
independent  model.  This  is  due  to  the  improved  modeling  of  the  scattering  physics.  Since 
the  current  required  for  the  extrapolation  is  only  computed  at  lower  frequencies,  large 
savings  in  computation  time  and  memory  can  be  achieved. 
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Fig.  1 .  Time-of-arrival  model  for  the  induced  current  at  point  P  accounts  for 

the  direct  incident  radiation  from  the  antenna  and  the  multiple  scattered 
waves  from  other  parts  of  the  platform. 


Fig.  2.  Extrapolation  of  a  simulated  signal  using  different  models. 


(a)  SNR  =  oo  (b)  SNR  =  lOdB 

Fig.  3.  Variation  of  Im(r„)  as  a  function  of  the  pre-multiplication  factor  cim  for 
the  simulated  signal  at  two  different  SNR  levels. 


Separation  between  two  arrival  times  (1/BW) 

Fig.  4(a)  RMS  error  in  estimating  the  time-of-arrival  tn  using  ESPRIT. 
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Separation  between  two  arrival  times  ( 1/BW) 
(b)  RMS  error  in  estimating  the  frequency  factor  a„. 


(c)  Percentage  RMS  error  of  the  extrapolated  signal  at  the  extrapolation  ratio  of  3:1. 

Fig.  4.  Error  performance  of  the  extrapolation  algorithm  as  a  function  of  the 
separation  between  two  arrival  times  at  different  SNR  levels. 
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Fig.  5.  2D  platform  geometry 
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Fig.  6(a)  Frequency  response  predicted  by  the  frequency-independent  time- 
of-arrival  model  at  d=  40°. 
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(b)  Frequency  response  predicted  by  the  frequency  dependent  model  at  0-  40' 


t  (ns) 


(c)  Time  domain  response  predicted  by  the  frequency  dependent  model. 
Fig.  6.  Frequency  extrapolation  example  for  the  2D  platform. 
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Fig  9.  (a)  Frequency  response  predicted  by  the  frequency  dependent 
model  at  <4/  =  30°,  <paz  =  -60°. 
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(a)  Computed  by  FISC 


(b)  Computed  through  extrapolation 


Fig.  10.  Comparison  of  the  radiated  field  generated  from  brute  force  FISC  computation 
and  frequency  extrapolation  as  a  function  of  frequency  and  azimuth  angle  at 
the  elevation  angle  of  <t>e\  -  50°. 
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ABSTRACT —  An  approximate-inverse  preconditioner  based  on  the  pre-defined  wavelet 
packet  (PWP)  basis  is  proposed  for  the  fast  iterative  solution  of  electromagnetic  integral 
equations.  The  PWP  basis  is  designed  to  achieve  a  sparse  representation  of  the  moment 
matrix  and  the  preconditioner  is  constructed  by  inverting  the  block-diagonal 
approximation  of  the  PWP-based  moment  matrix  and  transforming  the  results  into  the 
space  domain.  Numerical  results  show  that  the  PWP  preconditioner  is  effective  in 
accelerating  the  convergence  rate  of  iterative  solution  to  moment  equations.  It  is  also 
demonstrated  that  by  properly  designing  the  block-diagonal  matrix  and  computing  the 
matrix  elements,  the  total  computational  complexity  and  memory  costs  for  the 
preconditioner  can  be  kept  to  0(NlogN),  making  it  applicable  for  the  multilevel  fast 
multipole  method. 

Index  Terms  —  Computational  Electromagnetics,  Preconditioning,  Wavelet  Packet  Basis. 

I.  Introduction 

Iterative  algorithms  are  commonly  used  to  solve  large-scale  moment  equations 
resulting  from  electromagnetic  integral  equations.  The  computational  cost  of  iterative 
solution  is  proportional  to  both  the  moment  matrix-vector  multiplication  operation  and 


the  number  of  iterations  required  for  a  convergent  solution.  The  multilevel  fast  multipole 
method  (MLFMM)  has  been  developed  to  reduce  the  time  complexity  and  memory  of  the 
dense  matrix- vector  multiplication  operation  from  0(N2)  to  0(NlogN)  [1-3].  However, 
the  total  solution  time  remains  dependent  on  the  number  of  iterations  required  to  achieve 
an  accurate  solution.  If  the  scatterer  contains  near-resonant  structures  such  as  deep 
cavities,  the  moment  matrix  is  not  well  conditioned  and  the  convergence  rate  of  the 
iterative  solution  can  become  very  slow  [4,  5].  In  these  cases,  a  preconditioning  matrix 
[P]  is  desired  to  accelerate  the  convergence  rate.  We  consider  a  left-preconditioning 
system: 

[P][Z]  J  =  [P]  E  (1) 

where  [Z]  is  the  moment  matrix,  J  is  the  unknown  induced  current  vector  and  E  is  the 
excitation  vector.  To  achieve  the  best  preconditioning,  [P]  should  be  as  close  to  the 
inverse  of  [Z]  as  possible.  Hence,  [P]  is  called  an  approximate-inverse  preconditioner.  In 
this  work,  we  set  out  to  construct  an  approximate-inverse  preconditioner  [P]  that  satisfies 
the  following  conditions: 

(a)  [P]  is  an  accurate  approximation  to  [Z]'1  to  be  an  effective  preconditioner. 

(b)  The  total  computational  complexity  and  memory  requirement  to  construct  the 
preconditioner  are  of  0(NlogN). 

(c)  The  total  computational  complexity  and  memory  requirement  to  implement  the 
preconditioning  operation  in  (1)  are  of  O(NlogN). 

Conditions  (b)  and  (c)  are  set  to  keep  the  complexity  of  the  preconditioning  on  par  with 
that  of  the  MLFMM,  making  our  preconditioner  applicable  for  the  MLFMM.  The 
traditional  preconditioning  methods,  such  as  incomplete  LU  factorization  (ILU), 


Frobenius-norm-based  approximate  inverse  and  polynomial  preconditioners  could  be 
effective  [6,  7].  But  they  are  computationally  too  complicated  to  meet  conditions  (b) 
and  (c).  In  [5],  an  approximate-inverse  preconditioner  was  constructed  by  inverting  a 
block-banded  form  of  the  original  moment  matrix.  Low  complexity  was  maintained  by 
choosing  a  bandwidth  that  is  independent  of  problem  size. 

Our  approach  to  the  problem  is  to  derive  an  effective  preconditioner  meeting 
conditions  (a)-(c)  based  on  the  wavelet  packet  basis  [8].  With  the  conventional 
subsectional  basis,  [Z]  in  (1)  is  dense,  making  it  difficult  to  find  an  effective 
approximate-inverse  preconditioner.  However,  if  we  can  transform  the  moment  equation 
to  a  new  basis  to  make  the  transformed  moment  matrix  sparse  and  diagonal  dominant,  it 
will  be  much  easier  to  find  an  effective  preconditioner.  We  have  recently  proposed  the 
Pre-defined  Wavelet  Packet  (PWP)  basis  for  the  efficient  representation  of  moment 
matrices  [9,  10].  The  PWP  basis  was  designed  to  match  the  oscillatory  nature  of  the 
Green’s  kernel  in  the  integral  equation.  Numerical  results  showed  that  the  PWP- 
transformed  moment  matrices  are  strongly  diagonal  dominant,  and  have  about  OfNlogN) 
non-zero  elements  for  large-scale  problems.  In  this  paper,  we  shall  demonstrate  that  an 
effective  approximate-inverse  preconditioner  in  the  space  domain  can  be  derived  from 
the  PWP  basis  with  no  more  than  O(NlogN)  computational  or  memory  cost.  It  is 
worthwhile  to  point  out  that  there  have  been  some  recent  works  on  using  the  conventional 
wavelet  transform  to  precondition  linear  systems  for  iterative  solutions  [11-14]. 
However,  the  requirements  for  computational  cost  and  memory  are  at  least  0(N2)  to 
implement  the  wavelet  transform  and  the  preconditioning  operation.  In  this  work,  we 


overcome  the  0(N2)  complexity  and  memory  bottlenecks  by  computing  the  PWP 
preconditioner  directly  from  the  PWP  basis  functions  [10]. 

This  paper  is  organized  as  follows.  In  Section  II,  we  give  a  brief  review  of  the 
PWP  basis  and  its  representation  of  the  moment  matrices.  In  Section  III,  we  describe  the 
design  and  construction  of  our  approximate-inverse  preconditioner  in  the  space  domain 
from  the  sparse  PWP-based  moment  matrix.  We  show  that  the  total  cost  for  the 
construction  and  preconditioning  operation  can  be  kept  under  0(NlogN).  Numerical 
results  are  presented  in  Section  IV.  Finally  we  draw  some  conclusions  in  Section  V. 


II.  PWP  Basis  and  Its  Representation  of  Moment  Equations 

Wavelet  packet  basis  is  a  set  of  orthonormal  functions  derived  from  the  shift  and 
dilation  of  a  basic  wavelet  function,  and  can  be  considered  as  a  generalization  of  the 
conventional  wavelet  basis  [15].  From  the  point  of  view  of  function  space 
decomposition,  the  wavelet  packet  basis  space  is  generated  from  the  decomposition  of 
both  the  scaling  function  space  and  the  corresponding  wavelet  function  space.  Let  us 
assume  that  \j/(x)  is  the  wavelet  packet  basis  function  with  the  finest  spatial  resolution 
available  for  signal  representation.  Using  the  “two-scale  equations,”  we  can  express  the 
wavelet  packet  basis  functions  in  the  next  scale  as  [16]: 

vl(n)  =  '5Lv(k)h{ln-k) 


^\(n)  =  ^\ir(k)g(2n-k) 
k 


(2) 


where  (h(k)}  and  (g(k)}  are  the  impulse  responses  of  two  quadrature  filters  H  (low-pass) 
and  G  (high-pass),  respectively.  The  functions  in  the  next  scale  become  coarser  in  spatial 
resolution  and  finer  in  spectral  resolution  after  the  filtering  and  down-sampling  in  (2). 


The  same  procedure  can  be  applied  recursively  to  the  outputs  of  (2)  into  subsequent 
scales.  Conversely,  the  decomposition  results  in  (2)  can  also  be  used  to  reconstruct  the 
original  sequence  by  using  another  pair  of  quadrature  filters  P  and  Q,  which  are  the 
mirror  filters  of  H  and  G,  respectively.  This  reconstruction  procedure  can  be  expressed 
as: 

V(x)  =  'Z  P(x  ~  2k)Vo  (*) + J  ?(x  -  2k)y\  (k)  (3) 

k  k 

where  (p(k)}  and  (q(k)}  are  the  impulse  responses  of  P  (low-pass)  and  Q  (high-pass), 
respectively.  The  functions  become  finer  in  spatial  resolution  and  coarser  in  spectral 
resolution  through  filtering  and  up-sampling  in  (3). 

A  complete  and  orthogonal  wavelet  packet  basis  can  be  generated  from  a  spectral 
decomposition  tree  that  starts  from  an  initial  mother  function  \\i(x)  with  the  finest  spatial 
resolution  by  using  recursive  two-channel  filtering  and  down-sampling  in  (2).  It  can  be 
shown  that  the  decomposed  functions  at  the  outermost  branches  of  the  tree  satisfy 
orthogonality  and  completeness  for  any  decomposition  trees,  and  thus  constitute  a 
wavelet  packet  basis  set  [17].  In  this  context,  the  conventional  wavelet  transform  (CWT) 
basis  can  be  considered  as  a  special  case  of  the  wavelet  packet  basis  when  the 
decomposition  tree  grows  along  only  the  lowest  frequency  branch.  However,  we  have 
shown  in  [18]  that  the  structure  of  the  CWT  tree  is  not  optimal  for  representing  the 
moment  matrix  from  electrodynamic  integral  equations.  Instead,  we  have  proposed  a 
special  wavelet  packet  tree  for  electrodynamic  problems  [9,  10].  This  class  of  wavelet 
packet  bases,  which  we  term  the  pre-defined  wavelet  packet  (PWP)  basis,  has  a  spectral 
decomposition  tree  that  zooms  in  along  the  free-space  wave  number  ko.  The  motivation 
for  this  tree  structure  is  to  ensure  that  the  basis  is  well-matched  to  the  oscillatory  nature 


of  the  Green’s  function  kernel  in  the  integral  equation.  Figs.  1(a)  and  1(b)  show 
respectively  the  conventional  wavelet  decomposition  tree  and  the  PWP  decomposition 
tree  for  N=32.  In  the  PWP  tree,  the  center  frequency  of  its  deepest  branch  is  as  close  as 
possible  to  ko.  Detailed  discussion  on  the  construction  of  the  PWP  basis  can  be  found  in 
[10]. 

Once  the  PWP  basis  is  defined,  it  is  straightforward  to  transform  the  original 
moment  equation  into  the  new  basis  representation  as  follows: 


[Z]J  =  E 

(4) 

where 

[Z]~  [M]t[Z][M] 

(5) 

J  =  [M]tJ 

(6) 

E=[M]TE 

(7) 

[M]t  denotes  the  transpose  of  [M],  and  is  the  PWP  transformation  matrix.  It  changes  the 
original  expansion/testing  functions  from  the  standard  subsectional  bases  into  the  PWP 
bases.  Note  that  [M]T  is  equal  to  [M]'1  since  [M]  is  an  orthonormal  transformation. 

[Z],  J,  and  E  are  the  moment  matrix,  induced  current,  and  excitation  vector 
represented  in  the  PWP  basis,  respectively.  Our  numerical  results  from  [10]  showed  that 
the  transformed  moment  matrix  under  the  PWP  basis  is  very  sparse.  In  particular  the 
number  of  above-threshold  elements  in  the  transformed  moment  matrix  grows  at  a  rate  of 
0(N  '  )  for  small-sized  problems,  and  approaches  O(NlogN)  for  large-scale  problems. 

Furthermore,  the  significant  elements  in  [Z]  are  located  mainly  along  the  diagonal  or 
near-diagonal  positions.  Unfortunately,  the  exact  locations  of  the  significant  elements  are 
difficult  to  predict.  Consequently,  the  complexity  of  the  algorithm  in  solving  the 
complete  moment  equation  is  still  hampered  by  an  N2  computational  bottleneck,  since 


every  element  of  the  matrix  must  be  computed  to  determine  whether  it  is  small  enough  to 
be  discarded.  However,  for  the  proposed  preconditioning  application  in  this  paper,  we 
shall  show  in  the  next  section  that  an  effective  approximate-inverse  preconditioner  can  be 
constructed  by  computing  only  those  elements  within  the  diagonal  blocks. 

III.  PWP-Based  Approximate-Inverse  Preconditioner  in  the  Space  Domain 

We  will  first  outline  our  approach  to  constructing  the  PWP-based  preconditioner. 
Since  the  PWP  transformed  moment  matrix  [  Z  ]  is  very  sparse  and  diagonal  dominant, 
we  approximate  [  Z  ]  by  a  block-diagonal  matrix  [  Z  bd]  that  consists  of  the  near-diagonal 
terms  of  [Z].  The  approximate-inverse  preconditioner  [P]  defined  in  (1)  is  then 

constructed  by  inverting  the  block-diagonal  matrix  [  Z  bd]  and  transforming  the  resulting 
matrix  back  to  the  space  domain: 

[P]  =  [M][Zb„r1[M]T  (8) 

Therefore,  the  preconditioned  moment  equation  becomes: 

[M]  [ Z m]-1  [M]t[Z]  J  =  [M][Z m]-1  [M]tE  (9) 

The  inverse  of  the  block  diagonal  matrix  [  Z  bd]  is  simply  the  inverse  of  its  individual 
diagonal  blocks,  and  the  inverted  matrix  [Zbd]"1  remains  block  diagonal.  By  properly 
choosing  the  block  sizes  of  [  Z  bd],  we  can  limit  the  computational  cost  of  the  inversion 
while  maintaining  the  effectiveness  of  the  preconditioner. 

The  actual  implementation  of  the  preconditioner  will  now  be  considered.  If  we 
solve  equation  (9)  using  an  iterative  algorithm  such  as  the  conjugate  gradient  method,  the 
complexity  of  each  iteration  is  proportional  to  that  of  the  product  between  the  matrix 
[M],  [  Z  bd]’1,  [M] T  or  [Z]  and  a  vector.  The  product  [Z]  and  a  vector  can  be  implemented 


using  the  MLFMM  algorithm  with  OfNlogN)  complexity  and  memory  requirement. 
Because  the  PWP  transformation  matrix  [M]T  has  about  OfNlogN)  non-zero  elements, 
the  multiplication  of  [M]T  or  [M]  with  a  vector  is  also  of  OfNlogN)  in  complexity. 
Therefore,  if  we  can  limit  both  the  number  of  non-zero  elements  in  [Zbd]’1  and  the 
complexity  of  constructing  [Zbd]'1  to  OfNlogN),  the  total  cost  per  iteration  in  (9)  will 
remain  at  OfNlogN). 

Our  design  of  the  block-diagonal  matrix  [  Z  bd]  is  shown  in  Fig.  2.  The  blocks  are 
designed  so  as  to  capture  the  most  significant  elements  in  the  PWP-transformed  moment 
matrix  [  Z  ].  The  block-diagonal  matrix  consists  of  the  following:  (i)  a  block  centered 
about  the  spectral  frequency  ko  with  block  size  Mi,  (ii)  a  set  of  diagonal  blocks  in  the 
remaining  upper-left  region  with  block  size  M2,  and  (iii)  a  set  of  diagonal  blocks  in  the 
lower-right  region  with  block  size  M3.  The  number  of  diagonal  blocks  with  size  M2  is 
[fN/2-  Mi)/  M2],  and  that  with  size  M3  is  [fN/2)/  M3].  The  block  from  (i)  is  usually  the 
densest  part  of  [  Z  ],  since  its  elements  correspond  to  the  interactions  between  PWP  bases 
with  the  longest  spatial  extent  and  spectral  frequency  close  to  ko.  They  tend  to  have  the 
strongest  interaction  with  each  other.  The  blocks  from  (iii),  on  the  other  hand,  are  nearly 
diagonal,  since  their  elements  correspond  to  the  interactions  between  PWP  bases  with  the 
highest  spectral  frequency  and  the  shortest  spatial  extent.  Therefore,  we  generally  choose 
block  sizes  with  M]  >  M2  >  M3.  Furthermore,  as  N  is  increased  we  keep  the  sizes  M2  and 
M3  fixed,  but  let  Mj  grow  slightly  with  problem  size  at  a  rate  of  N  to  capture  even  more 
of  the  dominant  elements  in  [Z].  This  growth  rate  is  chosen  to  ensure  that  the  cost  of 
inverting  the  block  remains  bounded  by  N,  as  the  cost  of  the  inversion  is  proportional  to 
the  cube  of  the  matrix  dimension.  The  inverse  of  [Zbd]  is  equivalent  to  the  inverse  of 


each  of  the  blocks  of  [Zbd]  in  Fig.  2.  Therefore,  under  our  construct,  the  total 
computational  complexity  to  invert  [  Z  bd]  is  proportional  to  N  and  the  resulting  number  of 
nonzero  elements  in  [  Z  bd] 1  is  also  proportional  to  N. 

Next,  we  discuss  the  computation  of  the  required  elements  in  [  Z  M].  Each  element 
Z  (m,  n)  of  the  matrix  [  Z  ]  can  be  directly  expressed  as  [  1 0] : 

Z(m,n)  =  XZ^0')G0‘,y>„O')  (10) 

*  J 

where  G(i,  j)  is  the  free  space  Green’s  function,  i  and  j  represent  the  indices  of  the 
discretization  grid  in  the  space  domain,  and  \j/k  is  the  k-th  PWP  basis  function.  It  can  be 
seen  from  (10)  that  the  computational  complexity  for  each  matrix  element  is  proportional 
to  the  product  of  the  spatial  supports  of  the  two  corresponding  basis  functions.  The 
support  of  a  PWP  basis  function  is  related  to  the  depth  of  the  branch  in  the  PWP 
decomposition  tree.  Let  us  assume  that  the  basis  functions  corresponding  to  branches  in 
the  first  stage  of  the  PWP  decomposition  tree  have  spatial  support  L  (i.e.,  the  order  of  the 
wavelet  filter).  The  complexity  for  computing  a  diagonal  (or  near-diagonal)  element  is 
then  L2.  The  support  of  a  PWP  basis  function  doubles  if  it  corresponds  to  the  branches 
on  the  next  deeper  stage  in  the  decomposition  tree  [10,  16],  Therefore,  the  spatial  support 
of  the  basis  functions  generated  from  the  tree  branches  at  stage  k  is  2k''L,  and  the 
complexity  to  compute  an  element  from  basis  functions  in  the  same  stage  is  (2k_lL)2. 
Since  the  maximum  possible  depth  of  the  wavelet  decomposition  tree  is  Kmax=log2N,  the 
cost  for  computing  an  element  using  bases  from  this  maximum  depth  would  then  be  N2. 
However,  as  we  shall  discuss  in  the  next  section,  the  best  preconditioning  performance  is 
often  achieved  using  a  tree  with  depth  much  less  than  the  maximum  possible  depth. 
Similar  observation  has  also  been  reported  in  [11].  Therefore,  if  we  impose  that  the 


maximum  depth  index  of  the  PWP  decomposition  trees  is  a  small  number  K,  the  number 

^  •  9  If  1  9 

of  operations  needed  to  compute  an  element  in  [Zm]  will  be  between  L  and  4  '  L  . 
With  0(N)  elements  in  [ZM],  the  total  complexity  to  compute  the  matrix  [ZM]  is  about 
0(N)  operations. 

To  summarize,  we  have  designed  a  preconditioner  based  on  the  block-diagonal 
approximation  of  the  PWP-transformed  moment  matrix.  The  steps  to  construct  and  carry 
out  the  preconditioning  operation  entail:  (i)  computing  the  elements  in  the  diagonal 
blocks  of  [Zbo]  one  term  at  a  time  using  (10),  (ii)  inverting  [Zbd]  one  block  at  a  time  to 
arrive  at  its  inverse  [ZM]'',  and  (iii)  carrying  out  the  preconditioning  operation  by  a 
successive  set  of  sparse  matrix-vector  multiplication  operations  in  (9).  It  is  shown  that 
the  computational  complexity  of  the  each  of  the  three  steps  can  be  limited  to  within 
O(NlogN)  operations,  thus  making  the  proposed  preconditioner  compatible  with  the 
MLFMM  in  solving  moment  equations. 

IV.  Numerical  Results 

Two  two-dimensional  targets,  an  open-ended  inlet  and  a  bent  structure  shown  in 
Fig.  3,  are  chosen  to  test  our  PWP  preconditioner.  With  a  large  depth-to-opening  ratio 
(l/m),  the  moment  matrices  from  the  scatterers  are  ill-conditioned,  and  the  convergence 
rate  is  slow  if  an  iterative  solver  is  used  to  solve  the  systems.  For  the  inlet  structure,  if  the 
thickness  parameter  t  is  equal  to  zero,  the  electric  field  integral  equation  (EFIE)  is  use  to 
construct  the  moment  equation.  Otherwise,  the  combined  field  integral  equation  (CFIE) 
is  used.  For  the  bent  structure,  EFIE  is  always  used  to  find  the  solution.  The 
discretization  of  the  structure  is  chosen  to  be  fixed  at  0.1  wavelength  while  the  electrical 


size  of  the  structure  is  varied.  In  all  cases  the  incident  plane  wave  is  E-polarized  and 
makes  an  angle  of  45°  with  respect  to  the  inlet  opening. 

Fig.  4  is  the  moment  matrix  [Z]  represented  by  the  PWP  basis  for  the  inlet 
structure  with  dimensions  l:m:t=  15:1:1  and  N=256.  It  is  shown  in  logarithmic  scale 
without  any  thresholding.  We  find  that  the  most  significant  elements  in  the  transformed 
matrix  match  the  patterns  we  chose  for  the  PWP  preconditioner  in  Fig.  2.  The  densest 
portion  of  the  matrix  is  located  in  the  upper  left  region,  with  the  strongest  elements  near 
the  ko  spectral  position.  The  wavelet  filter  used  to  generate  the  PWP  basis  functions  is 
Daubechies  filter  with  order  16  (with  7  vanishing  moments). 

We  use  both  the  conjugate-gradient  squared  (CGS)  method  and  the  biconjugate 
gradient  stabilized  (BiCGSTAB)  method.  Both  of  them  are  effective  iterative  solvers  for 
linear  equations  with  nonsymmetric  coefficient  matrix  [19].  CGS  is  used  in  most  of  the 
cases  in  this  work  since  it  converges  faster  than  BiCGSTAB.  However,  we  use  the 
BiCGSTAB  to  compare  the  detailed  convergence  rates  of  the  different  preconditioners,  as 
the  convergence  behavior  of  the  CGS  is  not  very  smooth.  In  the  application  of  the 
iterative  solvers,  the  initial  current  guess  is  Jo=0,  and  the  stopping  criterion  used  is  when 
the  relative  residual  error  ||rn||/||r0|  <  1  O'6,  where  r„  =  E-[Z]Jn  and  r0-  E. 

We  first  examine  how  the  iteration  number  changes  as  the  electrical  size  of  the 
scatterer  is  increased  in  the  following  three  cases: 

Case  1:  The  scatterer  is  the  inlet  structure  shown  in  Fig.  3(a)  with  /  =  0  and 
l-m— 15:1.  The  moment  equation  is  formulated  with  the  EFIE.  For  the  PWP 

preconditioner,  we  choose  M2=32,  M3=l,  and  Mi  grows  at  a  rate  of  \[n  starting  from  32 


atN=128. 


Case  2:  The  scatterer  is  the  inlet  structure  with  /:w:f=15:l:l.  The  moment 


equation  is  formulated  with  the  CFIE.  The  PWP  preconditioner  is  constructed  the  same 
way  as  that  in  Case  1 . 

Case  3:  The  scatterer  is  the  bent  structure  shown  in  Figure  3(b)  with  l\m- 12:1. 
The  EFIE  is  used,  and  the  parameters  for  the  PWP  preconditioner  are  chosen  as:  M2= 

M3=  32,  and  M]  grows  at  a  rate  of  \[n  starting  from  32  at  N=128. 

The  iteration  numbers  required  versus  problem  sizes  are  shown  in  Figs.  5(a)-(c) 
for  Cases  1-3,  respectively.  We  increase  the  problem  size  by  proportionally  increasing 
the  electrical  size  of  the  scatterer  while  keeping  the  discretization  interval  fixed  at  0.1 
wavelength.  The  results  for  the  moment  equation  without  any  preconditioning  and  that 
with  the  PWP  preconditioner  are  plotted.  For  comparison,  we  also  show  the  results  from 
a  simple  block-diagonal  preconditioner  in  the  space  domain.  It  is  constructed  by 
inverting  the  block-diagonal  form  of  the  original  space-domain  moment  matrix.  The 
block  sizes  are  uniform  and  the  total  number  of  non-zero  elements  is  kept  the  same  as 
that  used  in  the  PWP  preconditioner.  We  find  that  the  iteration  number  grows  very 
rapidly  with  the  problem  size  if  no  preconditioning  is  applied.  With  the  PWP 
preconditioners  applied,  the  iteration  numbers  are  significantly  reduced  and  they  increase 
much  more  slowly  with  problem  size.  Furthermore,  the  PWP  preconditioner  performs 
better  than  the  equivalent  space  domain  preconditioner.  This  is  because  the  PWP 
preconditioner  captures  more  of  the  energy  in  the  original  moment  matrix.  Note  that  the 
CFIE  results  shown  in  Fig.  5(b)  have  smaller  iteration  numbers  than  the  results  from  the 
EFIE  cases  shown  in  Figs.  5(a)  and  5(c).  However,  the  effectiveness  of  the  PWP 
preconditioner  is  still  apparent.  In  Fig.  5(c),  we  also  plot  the  results  using  the 


conventional  wavelet  transform  (CWT)  basis  instead  of  the  PWP  basis.  The  CWT 
preconditioner  is  constructed  in  exactly  the  same  way  as  the  PWP  preconditioner,  except 
the  approximate  matrix  is  formed  from  the  elements  of  the  CWT  transformed  moment 
matrix  rather  than  the  PWP  transformed  one.  We  observe  that  the  PWP  preconditioner 
outperforms  the  CWT  preconditioner  since  the  PWP  basis  can  more  efficiently  represent 
the  moment  matrix.  Although  not  shown,  similar  comparisons  between  the  PWP  and 
CWT  preconditioners  are  also  found  for  Cases  1  and  2. 

Figs.  6(a)-(c)  show  the  convergence  behaviors  for  a  fixed  problem  size  (N=512) 
for  Cases  1-3,  respectively.  The  BiCGSTAB  algorithm  is  used  as  the  iterative  solver.  In 
each  figure,  the  residual  error  is  plotted  as  a  function  of  the  number  of  iterations  for  the 
moment  equation  without  any  preconditioning,  with  the  space  domain  and/or  the  CWT 
preconditioner,  and  with  the  PWP  preconditioner.  We  observe  that,  with  the  PWP 
preconditioning,  the  relative  residue  decreases  the  fastest.  In  addition,  the  convergence 
behaviors  are  more  stable.  Similar  convergence  behaviors  are  found  for  larger  problem 
sizes. 

Next,  we  investigate  the  effect  of  the  depth  index  K  in  the  wavelet  tree  on  the 
resulting  PWP  preconditioner.  For  a  scattering  problem  with  a  problem  size  of  N,  the 
maximum  depth  index  for  the  PWP  decomposition  tree  is  (log2N).  However,  this  does 
not  necessarily  lead  to  the  best  preconditioning  performance.  Fig.  7  plots  the  required 
iteration  number  versus  the  depth  index  K  used  in  constructing  the  PWP  preconditioner 
for  N=1024.  The  results  show  that  the  curves  flatten  considerably  after  K=5,  implying 
that  there  is  not  much  to  be  gained  beyond  this  depth  in  the  wavelet  tree.  The  reason  for 
this  phenomenon  is  because  although  the  maximum-depth  PWP  tree  generates  the 


sparsest  transformed  moment  matrix  [10],  the  long  PWP  basis  functions  resulting  from 
that  tree  shift  some  of  the  energy  into  the  non-diagonal  blocks,  making  the 
preconditioning  results  worse.  For  the  problem  sizes  we  have  studied  from  128  to  4096, 
the  optimal  PWP  tree  depths  are  usually  between  3  and  6.  Since  the  computational  cost 
of  constructing  the  preconditioner  is  proportional  to  the  maximum  depth  of  the  PWP 
decomposition  tree,  we  should  not  go  beyond  this  depth  when  constructing  the  PWP 
preconditioner. 

Finally,  we  examine  the  computational  cost  of  the  PWP  preconditioner.  Fig.  8 
shows  the  CPU  running  time  required  to  generate  the  PWP  preconditioning  matrix  [  Z  bd] 
for  Case  1  versus  problem  size  N.  The  direct  element  computation  in  equation  (10)  is 
carried  out  and  the  PWP  tree  depth  K  is  increased  slightly  from  3  to  6  as  the  problem  size 
is  changed  from  128  to  4096.  Fig.  9  shows  the  CPU  running  time  needed  to  implement 
the  multiple  matrix-vector  products  required  for  the  carrying  out  the  PWP 
preconditioning  in  (9).  We  exclude  the  time  needed  for  the  first  moment  matrix  and 
vector  product  in  (9),  as  that  operation  can  be  done  via  the  MLFMM.  We  find  that  the 
CPU  running  time  grows  at  rate  close  to  O(NlogN)  in  the  both  cases.  This  is  consistent 
with  the  estimates  given  in  Section  III  and  meets  the  objectives  we  have  set  forth  for  the 
preconditioner. 


V.  Conclusions 


A  preconditioner  for  the  moment  equation  based  on  the  pre-defined  wavelet 
packet  basis  has  been  proposed  in  this  paper.  Due  to  the  vanishing  moments  of  wavelet 


basis  functions,  the  PWP-based  moment  matrix  is  sparse  and  diagonally  concentrated. 
Consequently,  an  approximate-inverse  preconditioner  can  be  more  easily  designed  and 
constructed  than  that  based  on  the  original  space-domain  moment  matrix.  The  PWP 
preconditioner  is  constructed  by  inverting  a  block-diagonal  form  of  the  PWP-based 
moment  matrix  and  transforming  the  resulting  matrix  back  into  the  space  domain.  It  has 
been  shown  that  the  complexity  for  the  construction  and  preconditioning  operation  can  be 
kept  under  0(NlogN)  in  both  computational  cost  and  memory  requirement.  Our 
numerical  results  showed  that  the  iteration  numbers  for  the  PWP-preconditioned  moment 
equations  are  significantly  smaller  and  grow  at  a  lower  rate  than  those  without 
preconditioning  or  preconditioned  using  a  space-domain  preconditioner.  In  addition,  the 
PWP  preconditioner  also  performed  better  than  the  equivalent  preconditioner  derived 
from  the  conventional  wavelet  basis.  Application  of  this  concept  to  three-dimensional  (3- 
D)  problems  is  currently  under  investigation,  with  the  goal  of  combining  the 
preconditioner  with  the  3-D  MLFMM  to  obtain  solutions  of  large-scale  problems. 

Finally,  although  we  have  chosen  the  block-diagonal  matrix  as  the  underlying 
structure  of  our  approximate-inverse  preconditioner,  we  have  found  that  a  block-banded 
matrix  [5]  or  a  thresheld  transformed  matrix  [1 1,  20]  gave  better  preconditioning  results 
than  the  block-diagonal  form.  However,  the  inverses  of  those  matrices  are  dense, 
resulting  in  0(N2)  computational  bottlenecks  for  implementing  the  matrix-vector  product 
in  the  preconditioning  step.  A  worthwhile  topic  is  to  find  a  more  effective  approximate 
matrix  in  the  PWP  basis  domain  than  the  block-diagonal  one,  while  preserving  the 
OfNlogN)  total  computational  cost  and  memory  requirement. 
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Figure  1.  (a)  Conventional  wavelet  decomposition  tree; 
(b)  Pre-defmed  Wavelet  Packet  (PWP)  decomposition  tree. 


Figure  2.  The  pattern  used  to  construct  [  Z bd]  from  [ Z  ]. 


(a)  (b) 

Figure  3.  The  test  scatterers:  (a)  an  open-ended  inlet,  and  (b)  a  bent  structure. 
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Figure  4.  The  moment  matrix  represented  by  PWP  basis,  [  Z  ], 
for  the  inlet  with  N=256  (in  logarithmic  scale) 


Figure  5.  Iteration  numbers  vs.  problem  sizes  for  solving  moment 
equations  with  different  preconditioning  methods  for  the  scattering 
structures  in  (a)  Case  1,  (b)  Case  2,  and  (c)  Case  3.  The  problem  sizes 
are  increased  by  proportionally  increasing  the  scatterer  sizes. 


(C) 

Figure  6.  Convergence  behaviors  of  the  preconditioned  systems 
versus  iteration  numbers  for  PWP  and  other  preconditioning  methods  in 
(a)  Case  1,  (b)  Case  2,  and  (c)  Case  3  with  N=512. 


Figure  7.  Effect  of  the  PWP  decomposition  level  on 
iteration  number  for  Cases  1-3  using  the  PWP  preconditioner 

with  N=1 024. 
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Figure  8.  CPU  running  time  vs.  problem  size  for  forming  the  PWP 
matrix  [  Z  M]  using  ( 1 0)  for  Case  1 . 


Figure  9.  CPU  running  time  vs.  problem  size  for  implementing  multiple 
matrix- vector  products  in  (9)  for  the  PWP  preconditioning  in  Case  1 . 
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I.  Introduction 

In  antenna  design  and  analysis,  the  mounting  platform  can  have  a  significant  effect  on  the 
antenna  radiation  characteristics.  However,  rigorous  solution  of  the  radiation  problem  over  a  complex 
platform  is  very  time  consuming,  and  the  computation  complexity  increases  dramatically  as  the 
frequency  increases.  In  this  paper,  we  present  a  model-based  frequency  extrapolation  technique  with 
which  the  radiated  field  over  a  broad  band  of  frequencies  can  be  obtained  using  the  rigorously 
computed  results  at  low  frequencies. 

Our  approach  entails  three  steps.  First,  the  induced  current  on  the  platform  surface  is  computed 
at  low  frequencies  using  the  method  of  moments  (MoM).  Second,  we  apply  the  time-of-arrival  model 
to  the  current  on  each  basis  element  on  the  surface.  The  model  coefficients  are  obtained  using  super- 
resolution  algorithm  ESPRIT  [1  ].  Finally,  the  induced  current  at  higher  frequencies  is  computed  using 
the  model  and  the  radiation  characteristics  are  calculated.  This  approach  is  similar  to  that  present  in  [2] 
and  [3]  for  radar  signature  extrapolation.  Our  results  show  that  when  the  frequencies  and  the 
discretization  of  the  platform  are  properly  chosen,  the  radiated  field  at  higher  frequencies  can  be 
extrapolated  with  only  moderate  computational  cost. 


II.  Extrapolation  Methodology 

As  the  first  step  of  the  extrapolation  process,  the  induced  current  on  the  target  is  computed  at  a 
small  set  of  points  at  low  frequencies.  Once  the  current  over  the  surface  is  computed  at  low 
frequencies,  we  apply  the  time-of-arrival  model  to  each  of  the  current  element.  In  this  model,  we 
assume  that  the  total  current  is  induced  by  different  scattering  mechanisms,  as  shown  in  Fig.  1 .  Each 
of  the  incident  mechanism  has  a  distinct  arrival  time,  so  that  the  current  can  be  written  as 

7(<u)  =  (1) 

«»! 

where  co  is  the  angular  frequency,  N  is  the  total  number  of  incident  mechanisms  and  t„  is  the  arrival 
time  of  the  n*  incident  mechanism. 

Since  the  incident  mechanisms  correspond  to  scattering  from  different  parts  of  the  target,  the 
maximum  difference  in  t„  is  related  to  the  size  of  the  target.  Thus  the  sampling  rate  in  the  frequency 
domain  should  be  high  enough  to  distinguish  these  time  events  from  all  parts  of  the  target.  Based  on 
this  consideration,  we  approximately  constraint  the  sampling  rate  in  frequency  to 
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(2) 


where  c  is  the  speed  of  light  and  D  is  the  maximum  dimension  of  the  target. 

The  model  coefficients  an  and  /„  are  obtained  using  the  superresolution  algorithm  ESPRIT, 
which  is  based  on  the  model  that  the  signal  consists  of  a  sum  of  exponential  and  additive  white  noise. 
Given  a  sequence  with  M  samples,  the  algorithm  can  estimate  the  number  of  exponential  N  and 
determines  the  amplitude  and  period  of  each  exponential  term.  The  basic  requirement  in  the  number  of 
samples  is  M  >  IN  +  1.  Once  the  model  parameters  are  found,  the  induced  current  at  higher 
frequencies  can  be  computed  using  (1).  The  radiated  field  is  then  easily  obtained  from  the  extrapolated 
current. 


III.  Results 

As  an  example,  we  consider  a  2-D  structure  as  shown  in  Fig.  2  (a).  We  are  interested  in  the 
radiation  pattern  over  a  frequency  band  from  0.15  to  0.45  GHz.  To  obtain  the  data  for  the 
extrapolation,  we  compute  the  induced  current  at  10  frequencies  from  0.15  to  0.24  GHz.  Then  we  use 
the  ESPRIT  algorithm  to  obtain  the  model  coefficients  for  each  current  element  and  compute  the 
radiated  field  based  on  the  model.  Fig.  2(a)  shows  the  radiated  field  of  a  horizontally  polarized  line 
source  as  a  function  of  frequency  at  an  elevation  angle  of  40°.  The  dashed  curve  is  obtained  from  the 
model-based  extrapolation  while  the  solid  curve  is  the  reference  brute-force  solution.  We  observe  that 
the  extrapolation  algorithm  correctly  predicts  the  peaks  and  null  positions  in  frequency,  indicating  a 
good  estimate  on  the  times-of-arrival.  Fig.  2(b)  is  the  time  domain  response  obtained  via  an  inverse 
Fourier  transform  of  the  frequency  data.  The  first  large  peak  corresponds  to  the  specular  scattering 
from  the  flat  surface  and  the  second  large  peak  is  due  to  the  scattering  from  the  step  region. 

We  notice  from  Fig.  2(a)  that  the  model-predicted  field  matches  well  with  the  reference  field  at 
the  first  10  frequencies,  but  drops  below  the  computed  field  as  frequency  goes  higher.  This  is  because 
the  intensity  of  the  field  radiated  by  the  line  source  is  proportional  to  the  square  root  of  the  frequency. 
We  can  overcome  this  by  compensating  this  effect  before  doing  the  extrapolation.  Thus  the  time-of- 
arrival  model  becomes 


where  4co  should  be  replaced  by  a  for  3-D  situations.  After  the  compensation,  the  extrapolation  result 
is  further  improved,  as  shown  by  the  frequency  and  time  responses  in  Figs.  3(a)  and  3(b). 

Finally  we  look  at  the  radiation  problem  of  the  3-D  platform  in  Fig.  4(a).  The  source  is  a 
horizontally  polarized  dipole.  The  solver  used  is  FISC  [4],  which  is  a  3D  MoM  code  based  on  the  fast 
muitipole  method.  Similar  to  the  2-D  case,  we  use  the  computed  current  at  10  frequencies  from  0. 15  to 
0.24  GHz  to  extrapolate  the  data  to  0.45  GHz.  The  extrapolated  radiation  field  as  a  function  of 
frequency  at  the  elevation  angle  of  40°  is  shown  in  Fig.  4(b)  as  the  dashed  curve.  The  reference  brute- 
force  solution  is  plotted  as  the  solid  curve.  The  major  features  of  the  radiated  field  are  well 
characterized  by  the  extrapolation  algorithm. 
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IV.  Conclusion  and  Discussion 


As  we  have  seen  from  the  results,  the  frequency  extrapolation  technique  is  an  efficient  way  of 
obtaining  the  radiation  pattern  over  a  broad  band  of  frequencies.  Computation  time  is  reduced 
dramatically  since  the  current  is  rigorously  solved  only  at  low  frequencies.  We  improved  the  result  by 
compensating  the  frequency  factor  of  the  source  in  the  time-of-arrival  model.  This  indicates  that  a 
wrong  frequency  dependence  in  the  model  may  result  in  errors  in  the  model  coefficients.  In  addition  to 
the  frequency  factor  of  the  source,  the  frequency  dependence  of  different  scattering  components  of 
each  current  element  could  be  different,  due  to  different  scattering  physics.  Further  incorporation  of 
these  effects  should  further  enhance  the  accuracy  of  the  extrapolation. 
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(a)  Frequency  response 


»■ 


Fig.  2.  Frequency  extrapolation  for  2-D  problem  (H-Pol). 
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(a)  Geometry  of  the  3-D  platform 
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(b)  Frequency  response 

Fig.  4.  Frequency  extrapolation  for  3-D  radiation.  -  FISC 

- FISC  at  10  point  +  extrapolation 
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L  Introduction 

The  numerical  characterization  of  antenna  radiation  in  the  presence  of  a  large, 
complex  platform  is  usually  very  time  consuming.  The  problem  is  further 
compounded  when  the  radiation  pattern  is  desired  over  a  broad  frequency  band. 
In  this  paper,  we  present  a  frequency  extrapolation  technique  for  the  antenna 
radiation  problem  to  obtain  the  radiated  field  over  a  wide  band  of  frequencies 
from  a  limited  set  of  frequency  calculations.  Our  approach  is  similar  to  our 
previous  work  on  radar  signature  extrapolation  [1].  Given  a  CAD  model  of  the 
platform,  the  induced  current  on  the  surface  is  first  computed  using  a 
computational  electromagnetics  (CEM)  simulator  at  a  set  of  low  frequencies.  A 
time-of-anival  model  is  then  applied  to  each  current  element  and  the  super¬ 
resolution  algorithm  ESPRIT  [2J  is  used  to  calculate  the  model  coefficients.  The 
frequency  response  of  the  current  is  extrapolated  based  on  this  model  and  the 
radiated  field  is  obtained  over  a  wide  band  of  frequencies.  The  CEM  simulator 
used  in  our  work  is  the  multi-level  fast  multipole  solver  FISC  (3]. 

As  a  follow-up  to  the  frequency  extrapolation  algorithm,  we  also  set  out  to 
extract  a  sparse  point-radiator  model  of  the  antenna  radiation  pattern  in  the 
presence  of  the  platform.  We  assume  the  total  radiation  over  some  extended 
frequency  and  aspect  window  can  be  approximated  by  the  radiation  from  a  set  of 
radiation  centers.  The  location  and  amplitude  of  the  radiation  centers  are 
determined  by  using  a  matching  pursuit  algorithm  [4].  To  speed  up  the 
parameterization  time,  we  estimate  the  location  of  the  radiation  centers  by 
utilizing  a  Fourier-based  ASAR  (Antenna  Synthetic  Aperture  Radar)  algorithm 
[5]  we  have  developed  previously.  The  resulting  sparse  point-radiation  center 
model  can  be  used  for  real-time  reconstruction  of  complex  radiation  patterns.  In 
addition,  it  can  be  used  to  pinpoint  cause-and -effect  in  antenna-platform 
interaction. 

IL  Frequency  Extrapolation 

«.  An  example  of  the  antenna-platform  problem  is  shown  in  Fig.  1.  The  antenna 
is  assumed  to  be  a  horizontal  dipole.  We  first  compute  the  induced  current  at  a 
number  of  low  frequency  points  using  FISC  based  on  the  CAD  model  of  the 
platform.  We  assume  that  the  total  current  is  induced  by  different  scattering 
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mechanisms  where  each  of  the  incident  mechanism  has  a  distinct  arrival  time  [1]. 
Also  taking  into  account  the  fact  that  the  amplitude  of  the  radiated  field  from  a 
dipole  source  is  proportional  to  the  frequency,  a  current  element  can  be  described 
by  the  following  model: 

J(a)  *  o)  £anc~  (1 ) 

nml 

where  co  is  the  angular  frequency,  N  is  the  total  number  of  incident  mechanisms, 
and  On  and  /„  are  the  amplitude  and  arrival  time  of  the  n*  incident  mechanism, 
respectively.  The  model  coefficients  <z„  and  t„  are  obtained  by  parameterizing  the 
computed  currents  using  ESPRIT,  which  is  based  on  the  model  that  the  signal 
consists  of  a  sum  of  exponential  and  additive  white  noise.  Once  the  model 
parameters  are  found,  the  induced  current  at  higher  frequencies  can  be  computed 
using  (1).  The  radiated  field  is  theft  .easily  obtained  from  the  extrapolated  current. 

In  the  example,  the  induced  current  is  computed  from  0.15  to  0.24  GHz  at  an 
interval  of  0.01  GHz.  The  current  is  then  extrapolated  and  the  radiated  field  is 
computed.  The  extrapolation  result  is  compared  with  the  brute-force  reference 
data,  which  is  also  computed  by  FISC  over  the  frequency  band  from  0.15  to  0.45 
GHz  at  6  «  50°  and  *  =  0°,  as  shown  in  Fig.  2(a).  The  corresponding  time  domain 
response  is  plotted  in  Fig.  2(b).  The  matching  between  the  two  is  good  in  both  the 
frequency  and  time  domain,  demonstrating  the  effectiveness  of  the  extrapolation 
algorithm. 

IIL  Model-Based  Parameterization 

Next,  we  set  out  to  find  a  sparse  model  to  represent  the  radiated  field  from  the 
antenna- platform  configuration.  We  assume  that  the  radiated  field  can  be 
approximated  by  the  radiation  of  a  set  of  point  radiators,  each  having  a  frequency- 
aspect  behavior  described  by: 

EUkQ#)**  Ae'^e^*0 ********* ***+***•>  (2) 

where  k  is  the  wave  number,  (x&  ya  zo )  is  the  location  of  the  radiation  center.  The 
origin  of  the  above  basis  is  illustrated  in  Fig.  3.  The  matching  pursuit  algorithm  is 
applied  to  extract  the  point  radiators  one  at  a  time.  To  speed  up  the  extraction 
timr  of  the  matching  pursuit  algorithm,  we  estimate  the  location  of  the  radiation 
centers  by  utilizing  the  Fourier-based  ASAR  algorithm.  We  inverse  Fourier 
transform  the  radiated  field  with  respect  to  the  frequency  and  angles  to  generate 
an  ASAR  image.  The  point  with  the  highest  intensity  is  first  located  in  the  ASAR 
image  and  its  amplitude  and  coordinates  serve  as  an  estimate  of  the  strongest 
radiation  center.  We  next  zoom  in  on  the  precise  location  of  the  radiation  center 
via  a  fine  search.  We  then  subtract  the  contribution  of  this  radiation  center  from 
the  total  radiated  field  and  iterate  the  process  until  the  energy  in  the  residual 
signal  has  reached  a  sufficiently  small  leveL 
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We  demonstrate  this  concept  by  extracting  the  radiation  centers  from  the 
frequency-extrapolated  FISC  data  from  0.6  to  0.9GHz.  The  angular  range  of 
observation  is  over  a  30-degree  window  in  both  azimuth  and  elevation.  The  first 
10  extracted  radiation  centers  are  shown  in  Fig.  4,  with  their  amplitudes 
represented  by  grayscale.  We  observe  that  the  dominant  platform  scattering 
comes  from  the  edge  diffraction  from  the  right  edge  of  the  platform  and  the  comer 
structure  formed  by  the  cylinder  and  plate.  Note  that  the  radiation  center  due  to 
the  edge  diffraction  is  slightly  off  the  platform.  This  is  due  to  the  limited 
resolution  of  the  matching  pursuit  algorithm.  Once  the  sparse  representation  is 
generated,  the  radiated  field  can  be  rapidly  reconstructed  using  the  radiation 
centers.  The  original  radiated  field  at  **  =  0  is  plotted  in  Fig.  5(a)  as  a  function  of 
the  kj  and  kj,  where  k,  ky  and  kt  are  defined  in  (2).  The  field  reconstructed  from 
the  first  10  radiation  centers  is  shown  in  Fig.  5(b).  The  two  patterns  match  well 
over  both  frequency  and  angle. 

IV.  Conclusions 

We  have  applied  a  frequency  extrapolation  technique  to  the  antenna  radiation 
problem.  The  radiated  field  over  a  broad  band  of  frequencies- can  be  extrapolated 
efficiently  based  on  the  computation  result  at  a  limited  number  of  low  frequency 
points.  In  addition,  a  radiation  center  model  of  the  platform  radiation  can  be 
found  by  performing  the  matching  pursuit  algorithm,  liie  model  proves  effective 
as  the  radiated  field  can  be  accurately  reconstructed  with  very  little  computation 
effort. 
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Fig.  1.  Geometry  of  3-D  platform. 


Fig.  3.  Radiation  center  model. 


(a)  Frequency  domain  response  (dB) 


Fig.  2.  frequency  Extrapolation. 
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Fig.  4.  Extracted  radiation  centers. 


(a)  Original  radiated  field 
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(b)  Reconstructed  radiated  field 
Fig.  5.  Reconstruction  of  the  radiated 
field  from  the  radiation  center  model 
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I.  Introduction 

Recently,  a  number  of  iterative  solvers  based  on  fast  algorithms  have  emerged 
in  computational  electromagnetics.  These  solvers  have  much  higher 
computational  efficiency  than  traditional  approaches.  However,  for  wide-band, 
wide-angle  RCS  calculations,  the  solver  has  to  be  executed  repeatedly  for  each 
angle  or  frequency,  which  results  in  large  computational  expenses.  A  number  of 
interpolation  and  extrapolation  approaches  have  been  developed  to  generate  the 
RCS  curve  with  as  sparse  a  data  sampling  as  possible  [1J.  In  this  paper,  we 
address  the  interpolation  problem  in  frequency  and  angle  using  a  model-based 
approach.  Our  approach  is  the  adaptive  feature  extraction  (AFE)  algorithm.  It  has 
been  used  recently  by  us  to  eliminate  the  aliasing  effect  and  construct  ISAR 
image  from  unevenly  undersampled  measurement  data  [2].  Unlike  standard 
interpolation  algorithms  which  suffers  from  theNyquist  sampling  limitation,  AFE 
can  overcome  the  Nyquist  sampling  criterion  by  using  uneven  sampling.  The 
essential  idea  of  the  adaptive  algorithm  is  to  search  and  extract  out  individual 
scattering  features  one  at  a  time  in  an  iterative  fashion.  When  applied  to  the 
present  problem,  the  interference  between  different  scattering  features  can  be 
avoided.  After  all  die  main  features  are  extracted,  the  current  and  the  RCS  on  a 
denser  grid  of  sampling  can  be  interpolated  by  summing  the  contributions  from  all 
die  extracted  scattering  features.  The  AFE  algorithm  is  tested  using  numerical 
examples  for  both  1-D  frequency  interpolation  and  2-D  frequency-aspect 
interpolation. 


U.  1-D  Frequency  Interpolation  Algorithm 
In  simple  curve  fitting,  if  we  assume  that  the  far-field  frequency  response  of  a 
target  is  R(f)  and  the  field  at  sampling  frequencies  fis  known,  the  interpolated 
results  will  have  the  following  form, 

(i) 

i«i 

where  N  is  the  number  of  the  original  sampling  points  and  W(f)  is  the 
interpolation  basis  functions.  The  interpolation  will  have  good  performance  when 
the  original  function  is  a  slowly  varying  function  of  frequency.  According  to  the 
Nyquist  sampling  criterion,  the  sampling  interval  in  the  frequency  domain  should 
be  less  than  l/(2r„ux)  to  ensure  good  interpolation  results,  where  is  the 
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longest  time  delay  between  the  different  scattering  mechanisms  caused  by  the 
structure.  The  curve  fitting  procedure  can  also  be  carried  out  for  the  induced 
current  at  each  point  on  the  target  This  has  been  found  to  have  an  advantage  over 
the  field-based  interpolation  [3],  since  the  phase-demodulated  current  has  much 
looser  Nyquist  sampling  criterion.  As  a  result,  the  number  of  sampling  points 
needed  to  carry  out  the  interpolation  can  be  much  less  than  that  usiilg  a  fitting 
scheme  for  the  far  field.  Howev#*,  this  observation  is  not  always  tide  and  only 
works  for  targets  where  the  scattering  mechanism  is  dominated  by  physical  optics. 
In  this  work,  we  apply  the  AJFE  algorithm  to  the  cuncnt  to  deal  with  the  higher- 
order  interactions  in  more  general  targets.  Our  approach  is  based  on  the  following 
time-of-anival  model  of  the  induced  current 

J(/)='ZBftxp(-j2nftf)  (2) 

P 

This  model  gives  a  good  description  of  the  multiple  scattering  physics  at  high 
frequencies  [4,5].  The  number  of  terms  is  left  unspecified  and  will  be  determined 
by  the  algorithm.  To  determine  the  model  parameters  Bp  and  we  use  an 
iterative  procedure.  We  first  project  the  sampled  frequency  response  onto  the 
complex  conjugate  of  the  model  bases  in  (2): 

(3) 

N  M 

where  the  subscript  p  denotes  it  is  in  the  pth  stage  of  the  iterative  procedure.  It 
should  be  noted  that  the  sampling  f  is  intentionally  chosen  to  be  randomly 
distributed  in  the  frequency  band  to  avoid  the  ambiguity  in  deteimining  the 
strongest  feature.  Then  the  strongest  feature  in  this  stage  is  decided  by  exhaustive 
search  over  all  values  of  t, 

Bf  *  arg  max[Af,  (f)]  (4) 

Once  the  strongest  feature  is  found,  a  remainder  signal  is  produced  by  subtracting 
out  the  pth  feature: 

^,C//)  =  ^(/<)-^cxp(-y2^)  (5) 

After  all  the  amplitude  and  time-of-anival  parameters  are  known,  the  interpolated 
current  function  at  any  frequency  within  the  band  can  be  calculated  using  formula 
(2).  Therefore,  the  far  field  at  the  denser  frequency  sampling  can  be  obtained  by 
integrating  the  induced  current 

To  verify  the  above  analysis,  we  use  the  scattering  from  three  2-D  circular 
cylinders  as  an  example.  The  structure  is  shown  inFig.l(a).  The  original  RCS  is 
computed  at  71  points  from  0.3  GHz  to  0.65  GHz  using  the  method  of  moments. 
The  result  is  shown  as  the  solid  line  in  Fig.  1(b).  Next,  we  use  the  values  at  18 
equally  spaced  points  and  cany  out  a  current-based  interpolation  using  simple 
spline  fitting.  The  result  is  shown  as  the  dotted  line  in  Fig.  1(b).  Deviations 
between  the  two  results  can  be  seen.  While  the  current-domain  spline  gives  good 
prediction  of  the  physical  optics  current  on  the  target,  the  multiple  interactions 
mechanisms  are  not  well  interpolated.  Next,  we  use  the  AFE  algorithm  to  cany 
out  the  interpolation.  Instead  of  18  equal  spaced  points,  we  select  18  points 
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randomly  from  the  original  71  points.  The  result  is  plotted  as  the  dashed  line  in 
Fig.  1(b),  which  agrees  with  the  original  calculation  much  better  than  die  simple 
spline  interpolation.  This  comparison  is  further  shown  when  we  Fourier  transform 
the  frequency  responses  to  the  range  profiles  in  Fig.  1(c).  We  can  see  that  most  of 
the  features  from  the  AFE  interpolation  coincide  with  the  reference  result  while 
the  current-domain  spline  interpolation  gives  strong  artifacts  among  the  real 
scattering  features. 

ID.  2-D  Frequency-Aspect  Interpolation  Algorithm 
The  AFE  technique  can  be  extended  to  2-D  interpolation  in  both  frequency 
and  angular  domains.  We  replace  the  model  in  (2)  by  the  following  2-D 
scattering  center  model  [5]: 

J(Ovt'LBr  **P(-j—(x,  cos0  +  y,sin0  +  lp)]  (6) 
r  c 

where  lp  represents  the  time  delay  caused  by  the  higher  order  interactions.  This 
model  can  be  considered  as  an  extension  of  the  1-D  time-of-arrival  model  in  (2). 
To  show  the  payoff  of  the  2-D  interpolation  algorithm,  we  set  out  to  construct  an 
ISAR  image  of  a  cylinder-plate  structure  shown  in  Fig.  2(a).  Fig.  2(b)  shows  the 
reference  ISAR  image  generated  from  7 lx  81  sampling  points  in  the  frequency- 
aspect  plane.  The  grayscale  image  has  a  dynamic  range  of  40  dB.  To  construct 
the  image  using  the  AFE  algorithm,  we  choose  200  points  randomly  from  the 
original  5751  points.  The  ISAR  image  obtained  using  the  2-D  AFE  interpolation 
is  plotted  in  Fig.  2(c).  Most  of  the  features  of  the  interpolated  image  agree  well 
with  the  brute-force  calculation.  The  dynamic  range  of  the  interpolated  result  is 
limited  by  the  imperfection  of  the  2-D  model  in  (6)  and  the  errors  in  the  AFE 
procedure. 
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Fig.  1  (a)  Geometry  of  the  three-cylinder  target  Fig.2(a)  Geometry  of  the  cylinder-plate  target. 


Fig.  1(b)  RCS  vs.  frequency  plots  using  current  Fig.2(b)  ISAR  image  generated  using  brute 
domain  spline  and  AFE  interpolation.  force  MoM  calculation  at  575 1  points. 


Fig.  1(c)  Range  profiles  generated  from  the  Fig.2(c)  ISAR  image  generated  from 
interpolated  results  using  current  domain  spline  interpolated  result  using  AFE  with  200 
and  AFE.  calculated  points. 
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